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A robust  control  design  methodology  is  presented  for  systems  subject  to  ellipsoidal 
uncertainty.  Both  the  robust  analysis  problem  and  the  robust  synthesis  problems  are 
considered,  and  both  discrete  and  continuous  systems  are  considered.  For  continuous 
systems  subject  to  ellipsoidal  uncertainty,  a frequency  dependent  rational  function  is 
constructed  to  measure  the  stability  margin  at  each  frequency.  This  rational  function  may 
have  multiple  minima,  therefore  a frequency  search  must  be  carried  out  to  find  the  overall 
stability  margin  of  the  system. 

An  alternative  method  for  testing  robust  stability  is  developed  that  is  based  on 
constructing  a stable  transfer  function  whose  frequency  response  magnitude  is  equal  to  the 
stability  margin  of  the  system  under  consideration.  The  overall  stability  margin  is  then 
found  as  the  Hoo  norm  of  this  transfer  function.  The  construction  of  this  transfer  function 
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requires  performing  two  spectral  factorizations.  This  method  is  applicable  to  both  discrete 
and  continuous  systems. 

A robust  control  synthesis  procedure  for  systems  subject  to  ellipsoidal  uncertainty  is 
also  developed.  This  method  consists  of  constructing  a fixed  order  Youla  parameter  to 
achieve  a robustly  stable  closed  loop  for  the  specified  level  of  uncertainty.  The  appropriate 
coefficients  of  the  Youla  parameter  are  found  through  a constrained  quasi-con  vex 
optimization.  The  design  methodology  is  applicable  to  any  nominally  stabilizing  controller, 
but  the  specific  case  of  a predictive  controller  is  considered  in  this  work.  The  resulting 
controller  is  a predictive  controller  that  is  robust  with  respect  to  real  parameter  variations, 
and  it  also  retains  the  nominal  servo  performance  of  the  original  predictive  controller.  The 
design  procedure  is  easily  modified  to  allow  the  incorporation  of  integral  action  into  the 
robust  controller,  guaranteeing  the  offset-free  rejection  of  asymptotically  constant 
disturbances.  It  is  also  possible  to  incorporate  a slightly  restricted  form  of  unstructured 
uncertainty  into  the  design. 
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CHAPTER  1 
INTRODUCTION 


1.1  Motivation 

The  control  engineer  always  faces  the  problem  of  uncertainty  in  the  plant  model  used 
for  control  design.  This  uncertainty  arises  from  numerous  sources,  including  unmodeled 
or  neglected  dynamics,  use  of  lumped  parameter  models,  and  variation  of  model  parameter 
values.  To  compensate  for  the  fact  that  the  model  used  is  not  an  exact  representation  of  the 
actual  system  to  be  controlled,  it  is  common  practice  to  design  a controller  that  will  yield 
acceptable  performance  for  a set  or  family  of  plants  that  includes  the  specific  (nominal) 
model  used.  If  the  true  system  is  modeled  well  by  some  plant  in  this  family,  not 
necessarily  the  nominal  model,  then  the  controller  will  adequately  regulate  the  true  system. 
Robust  control  is  the  process  of  designing  a controller  for  a set  of  plants  rather  than  a single 
model.  This  design  technique  requires  the  solution  of  two  closely  related  problems:  the 
robust  analysis  problem  and  the  robust  synthesis  problem. 

The  problem  of  testing  whether  or  not  a closed-loop  system  remains  stable  when 
subject  to  perturbation,  called  the  robust  analysis  problem,  has  been  studied  since  the 
earliest  days  of  feedback  control  theory.  Indeed,  much  of  the  pioneering  work  of  this  field, 
such  as  that  done  by  Nyquist  (1932)  on  feedback  amplifiers,  considered  system  stability 
subject  to  parameter  variation.  The  classical  stability  margins-the  gain  margin  and  the 
phase  margin-are  consequences  of  this  concern  about  stability  robustness  of  closed-loop 
systems.  However,  the  only  available  analytic  tool  for  testing  the  stability  of  a 
characteristic  polynomial  was  the  Routh-Hurwitz  criterion.  The  stability  conditions  derived 
from  the  criterion  yielded  quite  unwieldy  algebraic  inequalities,  even  for  a small  number  of 
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uncertain  parameters,  so  it  appeared  the  task  of  deriving  succinct  robust  stability  conditions 
might  be  insurmountable.  Less  attention  was  given  to  this  problem  as  the  field  of  optimal 
control  matured.  The  linear  quadratic  regulator  (LQR)  design  developed  in  the  early  1960s 
seemed  to  be  capable  of  producing  controllers  that  "guaranteed"  good  stability  properties. 
However,  Doyle  (1978)  demonstrated  that  the  stability  "guarantee"  disappears  when  the 
LQR  is  implemented  as  output  feedback.  This  result  emphasized  not  only  the  importance 
of  solving  the  robust  analysis  problem,  but  also  of  solving  the  robust  synthesis  problem, 
that  is,  actually  designing  a controller  that  could  guarantee  stability  for  a family  of  plants. 

The  stability  analysis  of  systems  subject  to  real  parametric  uncertainty  began  to  receive 
renewed  attention  with  the  result  of  Kharitonov  (1979)  on  the  stability  of  interval 
polynomials.  This  important  result  states  that  the  stability  of  a general  characteristic 
polynomial  with  coefficients  that  vary  independently  within  real  intervals  can  be  checked  by 
testing  the  stability  of  only  four  specific  polynomials,  regardless  of  the  original  polynomial 
order.  Suddenly,  it  appeared  that  analytic  solutions  might  be  available  for  many  types  of 
real  parametric  uncertainty.  Interval  systems  received  the  lion's  share  of  the  resulting 
interest.  Several  strong  results,  such  as  the  sixteen  plant  theorem  of  Barmish  et  al.  (1992), 
and  the  generalized  Kharitonov  theorem  due  to  Chapellat  and  Bhattacharayya  (1989) 
followed,  but  these  restrict  the  allowable  controller  structure,  and  represent  only  analysis 
results.  A more  general  exact  method  for  analysis  of  real  parametric  uncertainty  is  given  by 
deGaston  and  Safonov  (1988),  however,  this  approach  requires  a large  numerical  effort  for 
computation  of  the  stability  margin.  It  is  noted  that  all  of  these  methods  assume  that  the 
allowable  parameter  variations  are  independent. 

For  systems  whose  parameters  are  interdependent,  fewer  results  are  available  for  the 
computation  of  the  stability  margin.  Soh  et  al.  (1985)  present  a method  for  calculating  the 
largest  stability  hypersphere  directly  in  the  coefficient  space.  This  technique  is  applicable  to 
continuous-time  or  discrete-time  characteristic  polynomials.  Biemacki  et  al.  (1987)  extend 
this  result  to  computing  the  stability  hypersphere  directly  in  the  plant  parameter  space,  for 
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multiple  input  or  multiple  output  systems.  Also  of  particular  significance  is  the  approach  of 
Tsypkin  and  Polyak  (1991)  who  provide  an  elegant  frequency  domain  method  for  testing 
the  robust  stability  of  a continuous-time  characteristic  polynomial  whose  coefficients  could 
vary  within  an  t?p-ball.  Similarly,  Guzzella  et  al.  (1991)  develop  a method  for  analyzing 
the  robust  stability  of  systems  subject  to  ellipsoidal  (weighted  l2)  uncertainty. 

The  ellipsoidal  uncertainty  description  is  a specific  class  of  parametric  uncertainty  that 
assumes  the  nominal  transfer  function  model  is  of  fixed  order  but  has  real  coefficients  that 
lie  inside  an  ellipsoid  in  the  coefficient  space.  This  uncertainty  description  is  ideally  suited 
to  the  case  where  parameter  estimation  techniques  are  used  to  determine  the  coefficients  of  a 
nominal  plant  transfer  function  and  also  a parameter  covariance  matrix.  The  matrix  can  be 
used  to  construct  an  ellipsoid  that  describes,  in  statistical  terms,  the  expected  values  of  the 
plant  transfer  function  coefficients.  The  works  of  Fogel  and  Huang  (1982),  Belfonte  and 
Bona  (1985),  and  Belfonte  et  al.  (1990)  consider  algorithms  where  parameter  estimates  are 
constrained  to  lie  inside  ellipsoidal  domains.  The  direct  incorporation  of  the  identification 
information  into  the  uncertainty  description  is  highly  desirable.  The  focus  of  this 
dissertation  is  the  robust  control  of  systems  subject  to  ellipsoidal  uncertainty. 

The  robust  synthesis  problem  was  first  studied  in  detail  for  the  case  of  unstructured 
uncertainty.  Here,  the  true  system  is  modeled  as  a nominal  transfer  function  matrix  subject 
to  additive  or  multiplicative  perturbation.  The  only  information  known  about  the 
perturbation  is  a norm  bound.  In  this  case,  the  robust  analysis  problem  is  straightforward 
to  solve.  For  unstructured  uncertainties  the  powerful  results  of  Hoc  theory  (Doyle  et  al. 
1989)  were  developed  to  solve  the  robust  synthesis  problem,  for  both  the  single-input, 
single-output  (SISO)  and  multi-input,  multi-output  (MIMO)  cases. 

In  most  situations,  more  information  than  simply  a norm  bound  is  known  about  the 
uncertainty  affecting  the  system;  this  is  referred  to  as  structured  uncertainty.  Ideally,  any 
available  information  about  the  uncertainty  should  be  incorporated  into  the  control  design  to 
reduce  conservatism.  This  added  information  often  makes  both  the  robust  analysis  and 
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synthesis  problems  difficult  to  solve  explicitly.  Perhaps  the  strongest  and  most  complete 
result  is  the  structured  singular  value  (p)  analysis  method  introduced  by  Doyle  (1982),  and 
the  almost  identical  multivariable  stability  margin  (1cm)  method  introduced  by  Safonov 
(1982).  For  any  linear  plant  subject  to  uncertainty,  the  p-analysis  method  involves 
constructing  a fictitious  nominal  system  transfer  matrix  and  a block  diagonal  matrix  whose 
blocks  correspond  to  the  actual  uncertainty  affecting  the  system.  Then,  the  value  of  p 
corresponds  to  the  smallest  uncertainty  that  will  destabilize  the  system.  This  value  is  also 
referred  to  as  the  stability  margin  of  the  system.  Unfortunately,  the  exact  computation  of  |l 
is  NP  hard  in  general  (Braatz  et  al.  1994),  and  usually  a convex  upper  bound  is  computed 
instead.  The  p synthesis  procedure  uses  the  upper  bound  to  compute  a controller  that  is  a 
local  minimizer  of  the  structured  singular  value.  Since  the  optimization  involved  is  not 
convex  in  both  of  the  variables  involved,  a local  optimum  is  the  best  that  can  be  guaranteed. 
The  presence  of  real  scalar  uncertainty  blocks  complicates  both  the  robust  analysis  and 
robust  synthesis  problems.  A convex  upper  bound  is  available  for  the  computation  of  p 
with  real  parameter  uncertainty  (Young  1994),  and  a controller  design  procedure  employing 
this  upper  bound  has  been  proposed  (Young  1996),  but  it  will  also,  in  general,  find  a local 
optimum.  The  generalized  structured  singular  value  (Chen  et  al.  1994a,  Chen  et  al. 
1994b)  is  a metric  similar  to  p but  allows  for  different  norms  to  be  used  on  the  real  and 
complex  blocks  of  the  uncertainty  matrix,  thereby  allowing  an  analysis  of  more  general 
uncertainty  descriptions.  Braatz  and  Crisalle  (1997)  have  extended  the  generalized 
structured  singular  value  to  include  ellipsoidal  uncertainty  descriptions.  However,  there  are 
no  controller  design  methodologies  based  on  the  generalized  structured  singular  value. 


1.2  Objectives 

The  first  goal  of  this  dissertation  is  the  development  of  a new  method  for  testing  the 
robust  stability  of  ellipsoidally  uncertain  systems.  A method  for  the  computation  of  the 
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robust  stability  margin  for  ellipsoidal  systems  is  developed  that  employs  a bisection  search 
instead  of  a frequency  search.  Two  spectral  factorizations  are  performed  to  construct  a 
stable,  real-rational  transfer  function  whose  magnitude  corresponds  to  the  system's  stability 
margin.  The  stability  margin  is  found  as  the  Hoo  norm  of  this  transfer  function,  allowing  a 
bisection  search  to  be  used  to  compute  the  stability  margin.  The  method  is  also  applicable 
to  both  continuous  time  and  discrete  time  systems. 

The  second  goal  is  to  develop  a robust  controller  synthesis  method  for  ellipsoidal 
systems.  A control  design  method  is  presented  that  is  applicable  to  any  ellipsoidally 
uncertain  system  with  a stabilizing  nominal  controller.  The  design  of  a robust  controller 
based  on  a nominal  predictive  controller  is  considered  specifically  in  this  dissertation. 
Predictive  control  is  a class  of  control  designs  that  uses  knowledge  of  the  future  set  point 
values  to  explicitly  predict  the  future  plant  output,  and  compute  a control  law  that  will  drive 
the  output  as  close  as  possible  to  the  desired  set  point.  Predictive  control  is  well- 
established  in  process  industries  (Seborg  1994)  mainly  due  to  its  flexibility  and  (relative) 
simplicity.  However,  few  results  are  available  for  the  robustness  analysis  of  predictive 
controllers,  especially  for  the  case  of  real  parametric  uncertainty.  Therefore,  the  design 
procedure  presented  here  addresses  a relevant  topic  in  robust  control. 


1.3  Structure  of  the  Dissertation 

The  dissertation  is  organized  as  follows.  In  Chapter  2 the  problem  of  robust  analysis 
of  ellipsoidally  uncertain  systems  is  introduced.  A necessary  and  sufficient  condition  for 
robust  stability  of  continuous  time  systems  subject  to  ellipsoidal  uncertainty  is  derived, 
using  frequency  domain  analysis.  The  condition  is  similar  to  that  proposed  for  discrete 
systems  by  Guzzella  et  al.  (1991).  The  problems  associated  with  this  method  are 
discussed,  and  a design  example  is  presented  to  illustrate  the  possibility  of  controller  design 
based  on  the  proposed  robust  stability  condition. 
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Chapter  3 presents  a new  method  for  stability  analysis  of  ellipsoidal  systems.  The 
method  consists  of  constructing  a stable  transfer  function  whose  frequency  response 
magnitude  is  equivalent  to  the  stability  margin  of  the  uncertain  system.  This  transfer 
function  is  constructed  from  knowledge  of  the  nominal  plant,  a nominally  stabilizing 
controller,  and  the  matrix  describing  the  uncertainty  ellipsoid.  Two  spectral  factorizations 
are  performed  to  determine  the  numerator  and  denominator  polynomials  of  the  transfer 
function.  The  Nyquist  robust  stability  margin  can  then  be  found  as  the  Hoc  norm  of  the 
transfer  function.  The  method  is  applicable  to  both  discrete  and  continuous  systems. 

In  Chapter  4,  the  construction  of  a controller  robust  with  respect  to  ellipsoidal 
uncertainty  is  developed.  The  underlying  method  is  due  to  Rantzer  and  Megretski  (1994), 
and  is  applicable  to  continuous  or  discrete  systems.  In  this  chapter,  the  discrete  case  is 
considered,  and  further,  the  nominal  controller  for  the  system  is  assumed  to  be  a predictive 
controller.  Robust  controllers  that  retain  the  nominal  performance  of  the  predictive 
controller  are  derived.  The  proposed  method  is  extended  so  that  the  resulting  robust 
controller  exhibits  integral  action,  assuring  offset-free  rejection  of  asymptotically  constant 
disturbances. 

Chapter  5 considers  a limiting-case  robust  control  synthesis  problem  for  ellipsoidal 
systems.  The  plant  is  assumed  to  be  first-order  and  the  controller  considered  is  a static 
gain.  This  case  is  of  interest  because  the  analysis  of  Chapter  2 is  not  strictly  applicable  and 
also  because  it  is  possible  to  directly  synthesize  a maximally  robust  controller  for  the 
system. 

Chapter  6 discusses  the  future  directions  of  this  dissertation.  Several  possible 
extensions  of  this  work  are  described,  along  with  a discussion  of  their  importance,  and 
possible  methods  for  solution  of  these  problems  are  given. 


CHAPTER  2 

STABILITY  ANALYSIS  FOR  ELLIPSOIDAL  SYSTEMS 


2.1  Introduction 

This  chapter  discusses  methods  for  stability  analysis  of  ellipsoidally  uncertain 
systems.  For  a continuous  time  system  subject  to  ellipsoidal  uncertainty,  a necessary  and 
sufficient  condition  for  robust  stability  is  derived.  For  discrete  time  systems,  the  necessary 
and  sufficient  condition  of  Guzzella  et  al.  (1991),  is  discussed  for  comparison.  In  both  the 
discrete  and  continuous  time  cases,  the  stability  condition  is  that  the  magnitude  of  a real- 
rational  function  of  frequency  be  less  than  one.  The  convex  nature  of  the  uncertainty 
regions  in  the  frequency  domain  allows  an  analytic  form  for  this  function  to  be  derived. 
This  function  can  have  multiple  local  maxima  and  minima;  a frequency  sweep  is  required  to 
find  the  global  maxima. 

A necessary  and  sufficient  condition  for  the  robust  stability  of  a continuous  time 
system  subject  to  ellipsoidal  uncertainty  is  given  in  Section  2.2.  The  results  of  this  section 
are  continuous  time  counterparts  of  the  results  given  in  Guzzella  et  al.  (1991),  although  an 
independent  derivation  is  provided.  Section  2.3  contains  a brief  review  of  the  discrete  time 
results  of  Guzzella  et  al.  (1991),  highlighting  the  differences  between  the  continuous  time 
and  discrete  time  cases.  Section  2.4  presents  an  example  that  uses  the  stability  margin 
defined  in  the  previous  sections  to  compute  a robust  controller  for  a water  heating  system, 
and  the  performance  of  this  controller  is  compared  to  controllers  designed  using  a nominal 
performance  measure. 
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2.2  Stability  Analysis  - Continuous  Time 

A derivation  of  the  robust  stability  analysis  problem  for  ellipsoidally  uncertain 
continuous  time  systems  is  given  in  this  section.  In  subsection  1,  the  general  ellipsoidal 
uncertainty  structure  is  introduced.  The  closed-loop  frequency  domain  uncertainty  regions 
that  arise  due  to  the  ellipsoidal  uncertainty  in  the  plant  are  discussed  in  subsection  2,  along 
with  the  derivation  of  a graphical  stability  test  for  ellipsoidally  uncertain  systems. 
Subsection  3 contains  the  necessary  and  sufficient  conditions  for  robust  stability  of  the 
continuous  time  ellipsoidal  uncertainty  problem.  Subsection  4 discusses  the  relationship 
between  the  parametric  stability  margin  and  the  Nyquist  robust  stability  margin. 


2.2.1  Ellipsoidal  Uncertainty 


The  general  plant  considered  in  this  analysis  is  the  linear,  strictly  proper  plant 

= fc„V'+...  + fcn  _ BW 

sn  + an_\Sn  1 + . . . + a0  A(.s) 
which  is  represented  by  the  coefficient  vector 

P = K-1  •••  a0  bn-\  ■■■bo]T  e 


(2.1) 


(2.2) 


The  values  of  the  coefficients  in  this  vector  are  uncertain,  however,  a nominal  value  of  the 


vector  is  assumed  known.  The  actual  value  of  the  parameter  vector  is  represented  by  the 
nominal  value  plus  an  additive  perturbation 

p = p°+5p  (2.3) 

This  perturbation  vector  is  constrained  to  lie  in  an  ellipsoid  Tp  in  the  parameter  space.  This 
ellipsoid  is  defined  by  a positive  definite,  symmetric  matrix  Qp  such  that 


%p  = {sP  E | 5ptq;'5p  < 1} 


(2.4) 


It  is  noted  that  the  numerator  and  denominator  polynomials  of  (2.1)  can  be  expressed  as  a 


nominal  polynomial  and  a perturbation  polynomial  as  follows 

, bV)  + AB(5) 

P{S)  = — rr 

Au(s)  + AAfO 


(2.5) 


where 
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AA(s)  = ' + . . . + 8oq 


and 


AB(s)  = $bn_[Sn  1 + . . . + 8b0 


The  coefficients  of  these  perturbation  polynomials  are  captured  in  the  vector  8p 


5p  = [5a„_,  ...  5a08Vi  S^0]T  e ^2" 


(2.6) 


An  example  of  an  ellipsoidal  uncertainty  region  for  a two  parameter  system  is  shown  in 
Figure  2.1. 


Figure  2.1.  Ellipsoidal  uncertainty  region  Figure  2.2.  General  feedback  interconnection 

Ellipsoidal  uncertainty  descriptions  arise  naturally  in  parameter  identification 
techniques  where  the  uncertainty  ellipsoid  is  associated  with  the  parameter-error  covariance 
matrix.  The  covariance  matrix  can  be  used  to  define  a matrix  Qp  that  corresponds  to  a 
certain  confidence  level  for  the  parameters.  Ellipsoidal  uncertainty  descriptions  are 
commonly  encountered  in  chemical  engineering  applications  where  model  parameters  are 
found  by  fitting  experimental  data  using  linear  or  nonlinear  regression  techniques.  Such 
uncertainty  descriptions  have  been  adopted  in  several  studies,  including  Biernacki  et  al. 
(1987)  where  ellipsoidal  domains  are  used  for  the  analysis  of  systems  characterized  with 
weighted  perturbation  bounds. 


+ 


v k 


Ellipsoidal  parametric  uncertainty  models  also  appear  in  various  other  contexts.  For 
example,  the  work  of  Agarwal  and  Bonvin  (1989)  on  Kalman  filtering  provides  a means 
for  estimating  the  covariance  matrix.  The  publications  by  Fogel  and  Huang  (1982), 
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Belfonte  and  Bona  (1985),  and  Belfonte  et  al.  (1990)  discuss  ellipsoidal  outer-bounding 
algorithms  which  produce  parameter  estimates  where  the  uncertainty  is  also  constrained  to 
lie  inside  ellipsoidal  domains  of  the  form  (2.4). 


2.2.2  Closed  Loop  Analysis 


Consider  the  problem  of  analyzing  the  robust  stability  of  a feedback  loop  containing  a 
fixed  controller  and  a plant  subject  to  ellipsoidal  uncertainty.  Figure  2.2  shows  the  general 
feedback  structure  adopted,  where  C(s ) is  a proper  controller  given  by 


C(s)  = 


^msn'  + ...  + P0 


P(s) 


(2.7) 


sm  +am_xsm  1 +...  + a0  a(s) 
and  is  represented  by  the  coefficient  vector  c 

c=K,_i ...  «0  P„, ...  Mt  s ‘*2”‘+1  <2-8> 

The  plant  P(s)  is  given  by  (2.1)  and  subject  to  the  ellipsoidal  uncertainty  described  by 
equations  (2.3)  and  (2.4).  It  is  assumed  that  the  plant  order  is  at  least  as  large  as  the 
controller  order,  i.e.,  n>m.  Furthermore,  it  is  also  assumed  that  the  index  k :=  n + m 
satisfies  k >2.  The  case  of  k = 1,  corresponding  to  a first  order  plant  ( n = 1)  and  a 
constant  controller  (m  = 0),  is  a limiting  case  that  is  discussed  in  detail  in  Chapter  5.  The 


characteristic  polynomial  of  the  feedback  loop 


G(5)  = G°(5)  + AG(5) 


(2.9) 


is  the  sum  of  nominal  and  perturbation  polynomials,  where 

G°(s)  = <x(s)A°(5)  + (3(s)B°(s)  = sk+  g°k_ + . . . + gS  (2.10) 

is  the  monic  nominal  characteristic  polynomial  and  the  perturbation  polynomial  is 

AG(5)  = gc(s)AA(j)  + P(  j)AB(s)  = + . . . + 5g0  (2.11) 

The  coefficients  of  the  characteristic  polynomial  define  the  characteristic  vector 

g = [a-!  -*o]T  S 9t*  (2.12) 

It  is  straightforward  to  show  that 

g = Sc  p + c 

where  Sc  is  the  Sylvester  matrix  for  the  controller  C(s),  and  has  the  form 
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sc  = 


1 0 
am-l  1 

a0 

0 

0 


0 


Pm 

0 Pm-1  P 


0 


1 Po 


a 


m- 1 


a0  0 


0 

0 

Pm 
Pm— 1 

Po 


e 91 


kxln 


and 


c = [oc m-\  •••  «o  0 •••  °]T  e * k 


The  nominal  value  of  the  characteristic  vector  is  generated  when  the  plant  parameter  vector 
assumes  its  nominal  value,  i.e.. 


g°  = scP° 


+ c = 


8k- 1 


8k- 2 


80 


(2.13) 


The  difference  between  the  actual  and  nominal  values  of  the  characteristic  vector  is  given  by 
the  perturbation  term 

5g:=g-g°  =Scp  + c-Scp°-c  = Sc(p-p°)  = Sc5p  (2.14) 

The  elements  of  8g  are  the  coefficients  of  the  polynomial  AG(s)  given  in  equation  (2.1 1). 
It  has  been  shown  in  Guzzella  et  al.  (1991)  that  if  5p  is  constrained  to  lie  in  the  set  Tp, 
then  8g  lies  in  the  ellipsoidal  uncertainty  region 


where 


5gTQg'5g<l} 


Qg  =ScQpSj  e XkKk 


(2.15) 

(2.16) 


is  a positive  definite,  symmetric  matrix. 

Throughout  the  remainder  of  this  chapter  it  is  assumed  that  the  controller  has  been 
chosen  so  that  it  stabilizes  the  nominal  system.  This  implies  that  G0(s)  corresponds  to  a 
Hurwitz  polynomial,  i.e.,  all  of  the  zeros  of  G°(s)  have  negative  real  part.  Therefore,  the 
problem  of  testing  if  a given  controller  robustly  stabilizes  the  system  becomes  the  problem 
of  testing  whether  any  allowable  AG(s)  produces  a non-Hurwitz  G(s).  This  problem  is 
analyzed  in  the  frequency  domain. 
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The  image  of  the  nominal  polynomial  G°(s)  as  s varies  along  the  imaginary  axis  is 

G°0'©)  = (sS  - g2 0)2  + •••)  + j(g\®  ~ S3®3  +•••):=  T*(e>)  + 7X/(co)  (2.17) 


Defining  the  vector 


it  follows  that 


where 


T(co)  := 


Tfl(CQ) 

LX/((0) 


G W 


T(CO)  = WT((0)g°  + t(C0) 


WT(co)  = 


... 

-to2 

0 l" 

wj(a>) 

... 

0 

to  0 

_wj((0)_ 

g * 


2\k 


(2.18) 


The  vector  t(0))  is  given  by 


and 


t(co)  = [(-l)*/2coA 
t(to)  = 


X 

o|  g ^K2  if  k is  even 


g 


if  k is  odd 


0 

The  image  of  the  polynomial  AG(s)  as  s varies  along  the  imaginary  axis  is  given  by 

AG(ja))  = (8g0  - 8g2co2  + ...)  + ;(8g,co  - 8g3co3  + ...):=  8x^(0) ) + j8x,((o)  (2. 1 9) 


and  can  also  be  expressed  as  the  vector 

Sx((0)  := 


SXfl(co) 

81,(0)) 


= WT(co)8g 


Guzzella  et  al.  (1991)  show  that  as  Sp  takes  on  all  possible  values  inside  the  parameter 
ellipsoid  Ep,  the  region  that  AG(/'co)  traces  out,  for  co  * 0,  is  an  ellipse  E^  in  the  complex 
plane.  This  ellipse  is  described  as 

?2 


where  the  matrix 


% w = |8x((o)  g 8xT(oa)Q[018T(to)  < lj 


Qco  = WT(co)Q„W(to)  g ST 


(2.20) 


(2.21) 


has  full  rank  for  all  nonzero  frequencies.  At  (0  = 0,  the  nominal  polynomial  is  simply 

G°(jO)  = go 


(2.22) 
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and  the  uncertainty  region  in  the  complex  plane  degenerates  into  a line  segment  defined  by 


(2.23) 


where  is  the  ( k,k ) element  of  the  matrix  Qg  defined  in  (2.16)  and  k = n + m. 
Uncertainty  regions  for  a non-zero  frequency  (0  and  for  to  = 0 are  shown  in  Figure  2.3. 

The  frequency  domain  uncertainty  set  is  called  a value  set,  or  template,  and  plays  a 
major  role  in  the  stability  analysis  of  ellipsoidal  systems.  This  analysis,  as  well  as  the 
classical  Nyquist  stability  criterion,  is  based  on  a result  from  complex  variable  theory  called 
the  principle  of  the  argument. 

The  principle  of  the  argument  states  that  as  a complex-valued  function  F(s ) is 
evaluated  along  a closed,  simple  curve  in  the  s-plane,  along  which  F(s)  has  no  zeros  or 

poles,  and  interior  to  which  F(s)  is  analytic,  then  the  net  change  in  angle  is 

Aarg{F(s)}  = 27t(np  -nz) 

where  nz  and  np  are  the  number  of  zeros  and  poles,  respectively,  of  F(s)  inside  the 
contour. 

The  principle  of  the  argument  can  be  used  to  test  the  stability  of  the  closed  loop  of 
Figure  2.2  by  applying  the  principle  to  the  characteristic  polynomial  G(s).  If  the  stability 
region  is  taken  to  be  the  open  left  half  of  the  complex  plane,  and  the  contour  is  the  classical 
Nyquist  contour,  traversed  in  the  clockwise  sense,  then  the  polynomial  G(s)  has  no  zeros 
in  the  right  half  plane  if  and  only  if 

Aarg{G(s)}  = 0 

since  a polynomial  has  no  finite  poles.  The  above  angle  change  is  the  net  change  as  the 
Nyquist  contour  is  traversed,  and  is  the  sum  due  to  the  semi-circular  part  and  the  imaginary 
axis  part.  The  change  in  angle  due  to  the  semicircular  portion  of  the  contour  is  -kn.  This 
follows  from  the  fact  that  on  this  part  of  the  contour  s = rejd  where  r — » and  0 varies 
from  y to— y.  Thus 

G(j)  ->  rkejke 
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and  the  net  change  in  angle  of  G(s ) is  -hi.  Therefore,  it  follows  that  the  change  in  angle  of 
G(/co)  as  to  varies  from  to  °°  must  be  kn  for  the  closed  loop  to  be  stable.  Since  the 
coefficients  of  G(s)  are  real,  the  plot  of  G(/co)  on  (-<*>, 0]  is  symmetric  to  the  plot  on  [0,-°°), 

fCK 

so  the  change  in  arg{G(/a>)}  on  [0,°°)  must  be  — . This  means  that  the  plot  of  G(/co) 

£ 

must  encircle  the  origin  - times  in  the  counterclockwise  direction.  This  important  result  is 
summarized  in  Lemma  2. 1 

Lemma  2.1  A polynomial  G(s)  of  degree  k has  all  of  its  roots  in  the  open  left  half 

k 

plane  if  and  only  if  the  plot  of  G(/co)  encircles  the  origin  of  the  complex  plane  — times  in 
the  counterclockwise  direction  as  to  varies  from  0 to 

A very  comprehensive  discussion  of  the  applications  of  the  principle  of  the  argument 
to  stability  testing,  including  a proof  of  the  Routh-Hurwitz  criteria,  is  given  in  Porter 
(1968). 

2.2.3  Robust  Stability  Analysis 

A robust  stability  test  can  be  constructed  from  the  result  of  Lemma  2.1.  The  nominal 
characteristic  polynomial  G°(s)  is  assumed  to  be  stable,  and  therefore  the  plot  of  G°(/co) 
has  the  correct  number  of  encirclements  for  stability,  as  described  by  Lemma  2.1.  Figure 

2.3  shows  the  uncertainty  ellipses  2^  at  several  frequencies.  The  band  that  is  swept  out  by 
these  ellipses  is  called  a Nyquist  envelope.  Each  point  in  a particular  ellipse  represents  an 
allowable  characteristic  polynomial  frequency  response  evaluated  at  that  particular 
frequency,  so  the  envelope  represents  the  frequency  response  plots  G(/'co)  for  all  allowable 
characteristic  polynomials.  If  the  envelope  does  not  include  the  origin,  then  all  the 
frequency  response  plots  have  the  same  number  of  encirclements  of  the  origin  as  the 
nominal  polynomial,  and  thus  all  allowable  characteristic  polynomials  are  stable. 
However,  if  the  Nyquist  envelope  contains  or  touches  the  origin,  then  at  least  one 
allowable  characteristic  polynomial  has  a different  number  of  encirclements  of  the  origin 
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than  the  nominal  polynomial  and  is  unstable.  Therefore,  the  system  is  robustly  stable  if 
and  only  if  the  value  sets  exclude  the  origin  at  all  frequencies. 

The  need  for  the  value  sets  to  exclude  zero  can  also  be  explained  in  terms  of  the  root 
locations  of  the  polynomial  G(s).  The  nominal  polynomial  G°(s)  has  all  its  roots  in  the  left 
half-plane  (LHP).  As  the  parameter  vector  p assumes  all  possible  values  in  Tp,  the  roots  of 
G(s)  move.  In  order  for  the  roots  to  travel  into  the  right  half-plane  (RHP),  they  must  cross 
the  imaginary  axis.  The  possibility  of  the  roots  moving  due  to  degree  dropping  has  been 
ruled  out  by  specifying  that  G(s)  is  monic.  Continuity  arguments  can  be  used  to  establish 
that  if  any  G(s)  has  roots  in  the  RHP,  then  an  allowable  G(s)  exists  that  has  at  least  one 
root  on  the  imaginary  axis.  However,  this  implies  that  G(/'co)  = 0 for  an  allowable  G(s) 
which  in  turn  means  that  the  origin  must  be  included  in  the  value  set  at  that  frequency. 
Therefore,  all  allowable  characteristic  polynomials  have  roots  in  the  LHP  if  and  only  if  the 
origin  is  excluded  from  every  value  set. 

Testing  this  zero  exclusion  condition  can  be  accomplished  using  the  critical  direction 
method  (Latchman  et  al.  1997),  modified  here  for  the  analysis  of  polynomials.  Consider 
the  following  quantities: 

i.  The  stability  segment  is  the  line  segment  joining  the  nominal  curve  G°(/(0)  and  the 
origin.  Note  that  the  length  of  the  stability  segment  is  simply  G0(yco)  . 

ii.  The  critical  direction  is  defined  as  the  direction  of  the  stability  segment,  and  is 
characterized  by  the  unit  vector 

G°0'a» 

d(;co)  = q 

G°0'a>) 

iii.  The  critical  perturbation  radius 

pc(co)  = maxja  e 91  | G°(;to)  + ad(y'co)  e 


These  quantities  are  illustrated  in  Figure  2.3.  It  is  noted  that  the  only  part  of  the  value  set 
important  for  robust  stability  is  that  part  that  lies  along  the  direction  to  the  critical  point. 
This  is  the  fundamental  premise  of  the  critical  direction  method. 
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Figure  2.3.  Stability  analysis  quantities 


As  stated  before,  the  closed  loop  is  unstable  if  any  value  set  includes  the  origin.  The 
value  sets  do  not  include  the  origin  provided  the  critical  perturbation  radius  is  smaller  than 
the  length  of  the  stability  segment  at  every  frequency.  Therefore,  a necessary  and 
sufficient  condition  for  robust  stability  is 


pc(co)  < G O'co) 


V®  g [0,°°) 


(2.24) 


Equation  (2.24)  is  particularly  useful  when  an  analytic  form  for  pc(CO)  is  available.  As 
Lemma  2.2  shows,  the  ellipsoidal  nature  of  the  value  sets  allows  such  a form  to  be  derived. 

Lemma  2.2  Suppose  the  nominal  characteristic  polynomial  G°(/a>)  given  by 
equation  (2.17)  is  subject  to  ellipsoidal  uncertainty  described  by  equations  (2.20)  and 
(2.23).  Then  the  critical  perturbation  radius  is 

||x(co)|L 

pc«0)=  , _2. <o>0  (2.25a) 

a/tT((0)Qco‘t(CO) 

Pc(0)  = V^U  0)  = 0 (2.25b) 


Proof:  The  proof  of  Lemma  2.2  is  given  in  the  Appendix. 
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By  definition  of  X(0),  and  from  equation  (2.22),  it  follows  that 

G°(;'g>)  =||x(co)||2 


G°(yO) 


So 


to  > 0 


m = 0 


(2.26a) 


(2.26b) 


The  results  of  Lemma  2.2  and  equation  (2.26)  can  be  used  to  write  the  necessary  and 
sufficient  conditions  for  robust  stability  given  in  (2.24)  as 


tT(co)  Qco'x(w) 


< ||<a>)|U 


go 


w > 0 


0 = 0 


or  they  can  be  cast  as 


xT(0)Q(1)1x(0) 


< 1 


go 


< 1 


0 > 0 


0 = 0 


The  Nyquist  robust  stability  margin  introduced  in  Latchman  et  al.  (1997)  for 
rational  transfer  matrices  can  be  reformulated  for  the  case  of  polynomial  systems.  Define 
the  frequency-dependent  Nyquist  robust  stability  margin  as 

1 


%(m)- 


'tT(0)Q(O1x(0) 


go 


0 > 0 


0 = 0 


(2.27) 


and  define  the  Nyquist  robust  stability  margin  as 

kN  = sup  %(m) 
oo>0 


(2.28) 


Theorem  2.1  Suppose  the  ellipsoidally  uncertain  plant  P(s)  described  by  (2.1)-(2.4) 
and  the  controller  C(s)  described  by  (2.7)  are  joined  in  unity  negative  feedback  as  shown  in 
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Figure  2.2,  with  the  nominal  system  stable.  Then  the  system  is  robustly  stable  if  and  only 
if  % < T 

Proof:  By  definition  of  both  Nyquist  robust  stability  margins  and  Lemma  2.2,  it  is 
immediately  apparent  that  the  necessary  and  sufficient  condition  given  in  (2.24)  is  satisfied 


Although  equation  (2.27)  gives  a closed-form  expression  that  can  be  used  to  test 
robust  stability,  the  frequency-dependent  Nyquist  robust  stability  margin  is  a rational 
function  of  frequency  that  may  have  local  minima.  This  implies  that  the  result  of  Theorem 


stability  margin.  While  each  point  of  the  search  requires  little  computational  effort,  it  is 
more  difficult  to  determine  the  range  of  frequencies  on  which  to  test  robust  stability.  The 
results  of  Chapter  3 provide  a method  for  testing  robust  stability  that  does  not  require  a 
frequency  sweep,  and  allows  the  value  of  t0  be  calculated  through  a bisection  search. 

2.2.4  Computation  of  the  Parametric  Robust  Stability  Margin 

The  parametric  stability  margin  is  defined  as  the  minimum  expansion  (or  contraction) 
of  the  uncertainty  region  in  the  parameter  space  required  to  bring  the  system  to  the  edge  of 
stability.  For  ellipsoidal  uncertainty,  the  parameter  space  uncertainty  region  is  contracted  or 
expanded  by  multiplying  the  matrix  by  a scalar  factor  a.  From  equation  (2.16),  it 
follows  that  multiplying  the  matrix  Qp  by  a scalar  a results  in  Qg  changing  to  aQg.  The 
new  uncertainty  region  in  the  characteristic  space  is  thus 


Using  this  definition  of  a scaled  uncertainty  set,  the  parametric  stability  margin  for 
ellipsoidally  uncertain  systems  can  be  expressed  as 


if  and  only  if  % < 1. 


V 


2. 1 must  be  tested  by  performing  a frequency  search  to  find  the  value  of  the  Nyquist  robust 


$ 

a - min  a(co) 
co 
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where 


a(to)  = min]  a | G(jco,5g)  = 0 
aeSR L 


5gea£^} 


The  relationship  between  the  parametric  stability  margin  and  the  frequency-dependent 
Nyquist  robust  stability  margin  is  detailed  in  the  following  lemma. 

Lemma  2.3  Suppose  the  ellipsoidally  uncertain  plant  P(s)  described  by  (2.1)-(2.4) 
and  the  controller  C(s)  described  by  (2.7)  are  joined  in  unity  negative  feedback  as  shown  in 
Figure  2.2,  with  the  nominal  system  stable.  Then  the  parametric  stability  margin  is 


* 

a 


Proof.  Let  the  frequency  that  defines  be  denoted  co*.  Then 

% =%((0*)  = 


tt(g)*)QC01*t(cd*) 


Multiplying  the  matrix  Qp  by  a scalar  results  in  Q?  changing  to  aQg  and  also  Qm  changing 


to  aQo  Thus,  the  system  with  the  scaled  uncertainty  satisfies 

1 Vot 


kN  =kN(co*)  = 


xT(a,*{l WUo)*)  ' '(“*) 


= fak 


■N 


Since  the  factor  a scales  all  the  values  of  the  frequency-dependent  Nyquist  robust  stability 
margin,  the  maximum  for  the  scaled  system  will  occur  at  the  same  frequency  as  that  for  the 


original  system.  The  parametric  stability  margin  is  defined  as  the  minimum  value  of  a that 
destabilizes  the  system,  which  corresponds  to  kN  = 1.  Therefore,  % = 1 = VoT*%  or 
a*  = kjj2 . V 


2.3  Stability  Analysis  - Discrete  Time 

The  stability  analysis  of  discrete  time  ellipsoidally  uncertain  systems  is  very  similar  to 
that  for  continuous  time  systems.  In  this  section,  the  few  differences  between  the  two 
cases  are  highlighted,  as  the  complete  derivation  can  be  found  in  Guzzella  et  al.  (1991). 
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In  the  discrete  time  formulation  of  the  ellipsoidal  stability  analysis  problem,  the  general 
linear  plant  considered  is  given  by 


b„_, + be  . B(z) 
Zn  + afl_izn~l  + ...  + OQ  ' A(z) 


(2.29) 


and  the  controller  is 


C(z)  = 


Pm*"1  + • • • + Po  P(^) 


(2.30) 


zm+am_lZm  ^+...  + ocq  a(z) 

The  plant  parameter  vector  p is  exactly  as  defined  in  (2.2)  and  the  uncertainty 
description  is  still  given  by  (2.3)  and  (2.4).  Also,  the  controller  vector  c is  given  by  (2.8). 
The  characteristic  polynomial  for  the  discrete  time  system  is  described  by  equations  (2.9)- 
(2.11)  with  the  Laplace  variable  5 replaced  by  the  discrete  time  variable  z.  The 
characteristic  polynomial  degree  remains  k:=n  + m.  Furthermore,  the  uncertainty  region 
“Eg  is  given  by  equation  (2.15)  and  the  uncertainty  matrix  Qg  is  given  by  (2.16). 

The  primary  difference  between  the  two  cases  is  the  frequency  domain  analysis  of  the 
characteristic  polynomial.  The  image  of  the  nominal  characteristic  polynomial  G°(z)  as  z 
varies  along  the  unit  circle  is  given  by 

G V“)  = (gS  + . . . + cos[(*  - 1)0 )]gf_!  + COSTCO]) 

+y(sin[a)]g1°  + . . . + sin  [(*  - l)co  + sin[£co])  (2.31) 

and  the  image  of  the  perturbation  polynomial  is 
AG(ejC0)  = (5g0  + . . . + cos [(k  - ) + jfsinfcoiS^j  + . . . + sin [(k  - l)co]5^A:_1 ) (2.32) 

These  polynomials  can  be  expressed  as  the  vectors 


T = 


Re 

G°(^t0)}' 

V 

Im  < 

G0(e;co)J 

xi_ 

= WTg°+t 


§x  = 


Re  {AG(e7'“)} 

"St/?" 

Im  {AG(ey(0)| 

8x7 

= WT5g 


and 
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where 

t = [cos(k©)  sin(kw)]T  6 9l2 


and 


WT  := 


cos[(k-l)©]  • 

• cos[©] 

f 

wj(©) 

sin[(k  - 1)©]  • 

• sin[©] 

0 

_w}(©)_ 

e 9* 


2 xk 


As  8p  takes  on  all  values  inside  "Ep,  the  uncertainty  region  AG( e JM)  that  is  generated 
around  G°( at  every  frequency  in  the  open  interval  (0,7t)  is  an  ellipse 


= jfo  e SR2 


5xTQC018x  < l} 


(2.33) 


where 


A(©) 

B(©) 

wJ(©)QgW/j(©) 

w^(w)Qgw/(©) 

_B(©) 

D(w). 

w]'(©)Qgw/?(©) 

w]'(©)Qgw/(©) 

(2.34) 


The  matrix  Qw  is  rank  deficient  at  the  two  frequencies  0 and  K.  At  these  points,  the 


uncertainty  regions  are  line  segments 


2b  = {8t*(0)  e SR  | |5t*(0)|  < ^wJ(0)Qgw*(0)  } (2.35) 

£*  = {&t*(tc)  e SR  | |5Tfl(7t)|  < ^Jwl(K)QgwR(K)  J (2.36) 

The  plot  of  all  possible  characteristic  polynomials  G(^)  for  co  e [7t,27t]  is  symmetric 
about  the  real  axis  to  the  plot  for  co  e [0,7i],  thus  only  the  range  to  e [0,7t]  need  be 
considered. 

The  principle  of  the  argument  is  used  to  deduce  the  root  locations  of  the  polynomial 
G°(z)  by  analyzing  the  number  of  encirclements  of  the  origin  of  the  plot  of  G°( as  © 
varies  between  0 and  k.  The  only  difference  from  the  continuous  time  case  is  that  now  the 
stability  region  is  the  interior  of  the  unit  circle  Izl  = 1.  The  counterpart  of  Lemma  2.1  for 
discrete  time  systems  is  stated  below. 
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Lemma  2.4  [Guzzella  et  al.  1991]  A polynomial  G(z)  of  order  k has  all  its  roots 
inside  the  unit  circle  if  and  only  if  the  number  of  counterclockwise  encirclements  of  the 


The  nominal  characteristic  polynomial  is  assumed  to  be  a Schur  polynomial,  i.e.,  it  has  all 
its  roots  in  the  unit  disc  Izlcl.  All  allowable  characteristic  polynomials  have  the  same 
degree  as  the  nominal,  so  the  plots  of  all  allowable  characteristic  polynomials  must  have  the 
same  number  of  encirclements  as  the  nominal  polynomial  for  the  loop  to  be  robustly  stable. 
Therefore,  for  robust  stability,  all  the  uncertainty  templates  must  exclude  the  origin.  As  in 
the  continuous  time  case,  the  highly  structured  nature  of  the  uncertainty  regions  allows  this 
zero  exclusion  requirement  to  be  tested  in  a direct  manner.  Figure  2.4  shows  the  plot  of  the 
nominal  polynomial  as  well  as  one  ellipsoidal  value  set  for  the  frequency  point  co*.  The 
value  sets  at  co=0  and  co=7t  are  line  segments.  The  stability  segment  at  the  frequency  co  is 
the  line  segment  joining  the  nominal  curve  G0(Vw*j  and  the  origin.  The  length  of  the 
stability  segment  is  simply  IG°fe/'a)jl.  The  critical  perturbation  radius  pc(C0)  is  the 


distance  from  the  center  of  the  value  set  to  its  boundary  along  the  stability  segment. 

V“*)  lm 


origin  performed  by  the  plot  of  G( e /c0)  as  co  varies  from  0 to  2k  is  equal  to  k. 


Figure  2.4.  Nominal  curve  G°(ei®),  value  sets  at  frequencies  co  = 0,  co*,  n,  the 
stability  segment  for  co*,  and  the  critical  perturbation  radius  pc(co*) 
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The  robust  stability  condition  for  the  discrete  time  case  is  based  on  the  knowledge  of  the 
explicit  form  for  the  critical  perturbation  radius,  which  is  given  in  Lemma  2.5 

Lemma  2.5  Suppose  the  nominal  characteristic  polynomial  G0^03)  given  by 
equation  (2.31)  is  subject  to  ellipsoidal  uncertainty  described  by  equations  (2.33),  (2.35) 
and  (2.36).  Then  the  critical  perturbation  radius  is 


Pc(03)  = - 

H®)ll2 

go  e (0,7t) 

(2.37a) 

Jtt(co)Qcd  c((o) 

Pc(0)  = 

tJ  wJ(0)Qgw^(0) 

co  = 0 

(2.37b) 

PcW  = 

^wJ(7t)QgW^(7t) 

© = n 

(2.37c) 

Proof:  The  proof  of  Lemma  2.5  is  entirely  analogous  to  the  proof  of  Lemma  2.2  and 
is  therefore  omitted. 


The  value  sets  % m do  not  include  the  origin  provided  the  critical  perturbation  radius  is 
smaller  than  the  length  of  the  stability  segment  at  every  frequency.  Therefore,  a necessary 
and  sufficient  condition  for  robust  stability  is 


Pc(C0)  < 


G°(ejU))  V to  e [0,7c] 


(2.38) 


Lemma  2.6  Suppose  the  conditions  of  Lemma  2.5  hold.  Then  the  condition  (2.38) 
is  equivalent  to  the  three  conditions 

<1  to  e (0,7t)  (2.39a) 


x (co)Qa)x((o) 


wJ?(0)Q„w/j(0) 


M<»l 


< 1 


co  = 0 


(2.39b) 


Jwl(n)Q  wR(n) 

; — < 1 

M*) 


to  = n 


(2.39c) 
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Proof:  The  form  of  pc(co)  is  given  in  Lemma  2.5.  For  0)  e (0,71), 


by  definition.  For  the  extreme  points  CD  = 0 and  co  = n,  T/(cd)  = 0,  and 
G°(e-/C0)  = |Tfl(co)|.  Plugging  these  values  into  (2.38)  yields  (2.39)  immediately.  V 

Defining  the  frequency-dependent  Nyquist  robust  stability  margin  as 

1 


*n(“)  = 


xT(co)Q<01x(ca) 


w«(0)QpW/j(0) 


M°)| 

|X^(7T)| 


co  e (0,7t) 


CD  = 0 


CD  = 7t 


(2.40) 


and  the  Nyquist  robust  stability  margin  as 

kN  = sup  &N(co) 

coe[0,7i] 

then  the  necessary  and  sufficient  condition  for  robust  stability  of  ellipsoidally  uncertain 
systems  is  given  in  Theorem  2.2. 


Theorem  2.2  Suppose  the  ellipsoidally  uncertain  plant  P(z ) described  by  (2.29), 
(2.2)-(2.4)  and  the  controller  C(z ) described  by  (2.30),  are  joined  in  unity  negative 
feedback,  with  the  nominal  system  stable.  Then  the  system  is  robustly  stable  if  and  only  if 

kN  < 1 (2.41) 


Proof.  Lemma  2.6  shows  that  the  satisfaction  of  the  three  conditions  (2.39)  is  the 
necessary  and  sufficient  condition  for  robust  stability.  However,  the  quantities  on  the  left 
hand  side  of  (2.39)  are  identical  to  those  on  the  right  hand  side  of  (2.40).  It  follows 
immediately  that  the  necessary  and  sufficient  condition  for  robust  stability  is  that 
/:n(cd)  < 1 for  CD  6 [0,;t],  or  < 1,  which  is  equation  (2.41).  V 
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2.4  Design  Example  for  Ellipsoidal  Systems 

In  this  section,  a simple  water-heating  control  system  with  a discrete  PI  controller  is 
used  for  illustrating  the  analysis  and  synthesis  methods  introduced  in  this  chapter.  The 
nominal  plant  parameters  are  used  to  tune  three  candidate  PI  controllers  using  standard 
techniques,  and  then  the  robustness  of  each  controller  is  analyzed.  It  is  verified  that 
controllers  cannot  be  robustified  from  knowledge  of  the  nominal  process  alone.  Finally,  a 
robust  control  design  based  on  the  results  of  Section  3 is  realized  via  numerical 
optimization. 


2.4.1  Plant  Model  and  Ellipsoidal  Uncertainty 


The  two-tank  system  with  recycle  shown  in  Figure  2.5  is  used  as  a model  for 
illustrating  the  robust  analysis  and  synthesis  techniques.  The  control  objective  is  to 
maintain  the  temperature  of  the  second  tank  (T2)  at  a desired  set  point  by  manipulating  the 
power  (R)  delivered  by  the  heater  located  in  the  first  tank.  The  actuation  on  the  heater  is 
performed  through  a zero  order  sample-and-hold  element.  The  only  available  measurement 
is  temperature  T2. 


Fi  F2 


1 r 

\ Tl 

Fi  +F2 

t2 

Fi+F2 

F, 

Fi  = 0.050  m^/min 
F2  = 0.150  m-Vmin 

p = 1000  kg/m^ 

Cp  = 4286  J/kg-C 
V = 1 


Figure  2.5.  Two  mixing  tanks  arranged  in  cascade  with  recycle  stream. 


The  actual  plant  model  is  generated  by  carrying  out  an  energy  balance  in  each  tank. 
Assuming  that  the  liquid  volume  V remains  constant  and  equal  in  both  tanks  one  arrives  at 
the  representation 


T2(s)  _ kj 
R(s)  x2s2  +2xs  + k2 


(2.42) 


where 
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r 

k,  = = 0.07  — 

(F,+F2)pCp  kW 

k2  = Fl  = 0.25 
F,+F2 

V 

x = = 5 min 

Fj  +F2 

The  actuation  on  R is  carried  out  by  means  of  a sample-and-hold  element  with  sampling 
period  T=100  sec.  The  corresponding  discrete  transfer  function  representation,  the  actual 
plant,  is  then 

T2(z)=  0.003  lz  + 0.0025 

R(z)  z 2 - 1.4932z  + 0.5134 

The  identification  of  the  nominal  plant  model  is  carried  out  using  a standard  least- 
squares  parameter  identification  method.  The  input-output  data  needed  for  identifying  the 
nominal  plant  G(z;  p°)  and  the  uncertainty  model  “Ep  is  generated  by  applying  a pseudo- 
random binary  input  to  the  heater  and  recording  the  temperature  response. 

In  the  simulation  it  is  assumed  that  the  temperature  data  are  acquired  by  a thermocouple 
which  perturbs  the  measurements  with  an  additive  white  Gaussian  noise  with  variance 
g2=4.5  10-3  C2.  A total  of  125  samples  are  gathered  using  a sampling  period  T=100  sec. 
The  collected  data  is  then  regressed  using  the  ARX  model 

T2  U)=  b'Z  + b°  R(z)  + e(z) 
z + ajz  + a0 

where  e(z)  is  assumed  to  be  a white-noise  sequence  of  unknown  variance.  Hence,  the 
plant  model  and  parameter  vector  for  this  representation  are,  respectively, 

P(z; P)=  b'Z  + bQ — , n = 2 (2.44) 

z + aiz  + a0 
and 

T 

P = h a0  bi  Z?0] 

The  least-squares  procedure  minimizes  the  functional  J(p)  = (y-Hp)T(y-Hp),  where 
y is  a vector  of  temperature  measurements  and  H is  the  regression  matrix  (Draper  and 
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Smith  1981).  The  minimum  is  realized  by  the  solution  vector  For  the 

two-tank  plant  under  study  the  solution  is 

p°  =[-1.1871  0.2087  0.0028  0.0038]1 
yielding  the  identified  nominal  plant  model 


0x  0.0028z  + 0.0038 
P(z]  p ) = — 


(2.45) 


- 1.1871z  + 0.2087 

An  ellipsoidal  uncertainty  model  is  readily  available  as  a by-product  of  the  least- 
squares  identification  technique.  A standard  result  (Draper  and  Smith  1981;  Crisalle  and 
Bonvin  1991)  states  that  a 100(l-a)%  joint  confidence  region  for  the  parameter  estimate  p° 
is  given  by  the  ellipsoidal  domain  given  in  (2.4) 

>2n 


with 


= {5peS\2''  5pTQ/?15p<l} 


Q„  = 2/iFEs' 


(2.46) 


In  the  above  equations,  2 n is  the  dimension  of  parameter  vector  p°,  k is  the  number  of 
sample  points  used  in  the  identification,  F = F(2n,  k-2n;  a)  is  a factor  characterizing  the 
1-a  quantile  of  an  J distribution  of  order  (2n,k-2n),  s2=(t-Hp°)T(t-Hp°)/(k-2n)  is  an 
estimate  of  the  noise  variance,  and  Zs  2=(HtH)  -1s2  is  an  estimate  of  the  covariance  matrix 
for  the  parameter  error.  The  statistical  interpretation  of  the  parameter  uncertainty  region 
provides  a sound  theoretical  justification  for  the  ellipsoidal  uncertainty  description. 

An  ellipsoidal  uncertainty  model  for  the  two-tank  process  is  obtained  by  substituting 
n = 2,  along  with  the  result  s2=  1.14  10'2C2  and  the  numerical  values  for  matrix  Z,  in 
(2.46)  to  get 


Q„  = F 


1.26xl0-2 

-1.23xl0-2 

-l.OOxlO-5 

3.51xl0-5 

-1.23xl0-2 

1.27xl0-2 

8.51xl0-6 

-3.62xl0-5 

-l.OOxlO-5 

8.51xl0-6 

4.25xl0-7 

-3.41xl0-8 

3.51xl0-5 

-3.62xl0-5 

-3.41xl0-8 

5.17xlCT7 

(2.47) 
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Hence,  using  the  above  form  for  the  matrix  in  (2.4)  represents  an  ellipsoidal  uncertainty 
model  which  can  be  interpreted  as  a region  with  a confidence  level  of  100(l-a)% 
depending  on  the  value  of  the  F-factor  chosen. 

The  fundamental  uncertainty  information  is  in  fact  contained  in  matrix  £ whose 
eigenvectors  determine  the  principal  directions  of  the  ellipsoid  Tp.  The  remaining  scalar 
factors  adjust  the  volume  of  the  ellipsoid  without  perturbing  the  principal  directions. 
Setting  the  F factor  equal  to  2.37  leads  to  a 95%  confidence  region;  however,  it  is  also 
possible  to  choose  the  value  of  the  F factor  to  encompass  a family  of  uncertainties  of  larger 
or  smaller  scope.  We  adopt  the  value  F=10  for  our  design,  thus  requiring  that  the 
controller  be  able  to  stabilize  a family  of  plants  belonging  to  a relatively  large  ellipsoidal 
region. 


2.4.2  Controller  Design 


The  controller  is  a discrete  proportional-integral  (PI)  compensator 

C(z;c)=PlZ  + f°  , c = [-1  p,  P0]T,  m = 1 

z - 1 

whose  associated  Sylvester  matrix  Sc  and  vector  c are 


0 

1 

-1 


Pi 

Po 

0 


0 

Pi 

Po. 


and  c = [-1  0 0]T 


For  the  nominal  plant  and  uncertainty  description  just  described,  three  PI  controllers 
are  synthesized  using  a conventional  tuning  method.  The  designs  proposed,  denoted  Cl, 
C2,  and  C3,  are  characterized  by  the  tuning  parameters  P]  and  P0  shown  in  Table  2.1.  The 
control  parameters  for  design  Cl  are  found  using  the  estimated  model  (2.45)  and  the  tuning 
settings  pi=Kc(l+0.5T/Ti)  and  Po=-Kc>  where  Kc  and  Ti  are  the  control  gain  and  integral- 
mode constants  of  a corresponding  analog  PI  controller  (Seborg  et  al.  1989).  In  order  to 
find  appropriate  constants  for  the  analog  controller,  the  step-response  of  (2.42)  is  first 
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approximated  by  a first-order  lag  and  a pure  delay  equal  to  one  half  of  the  sampling  period. 
Parameters  Kc  and  Tj  are  then  determined  from  the  correlations  minimizing  the  ITAE 
(Seborg  et  al.  1989). 

Designs  C2  and  C3  are  intended  refinements  to  C 1 based  on  the  observation  that  the 
[3 1 parameter  for  the  latter  is  near  a nominal  stability  limit.  In  fact,  when  Po=-98,  a 
nominally  stable  loop  is  obtained  for  98<(3i<149.  Thus,  the  higher  values  of  parameter  Pi 
in  designs  C2  and  C3  represent  an  attempt  to  move  the  controller  parameters  away  from  the 
nominal  stability  boundary.  Clearly,  the  proposed  tuning  refinements  are  made  based  on 
only  nominal  stability  considerations.  Figure  2.6a  shows  that  all  three  controllers  produce 
nominally  stable  closed  loops. 


Figure  2.6.  Closed-loop  responses  of  conventional  controllers  Cl,  C2,  and  C3, 
and  of  the  optimally  robust  controller  C*  to  a 5 °C  step  change  in  set- 
point.  (a)  Nominal  response,  (b)  actual  response. 


The  stability  robustness  of  all  three  candidate  controllers  is  then  analyzed  using  the 
critical  direction  technique.  Values  for  the  robustness  parameter  % are  calculated  by  means 
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of  an  exhaustive  numerical  search  using  a frequency  increment  Aco=0.01.  The  result 
obtained  is  given  in  the  last  column  of  Table  2.1.  It  is  concluded  that  the  loop  involving  C 1 
is  robustly  stable  (&n<1)>  but  that  neither  controller  C2  or  C3  is  robustly  stable.  In 
particular,  controller  C3  produces  an  unstable  loop  when  it  is  used  to  control  the  actual 
plant,  as  shown  in  Figure  2.6b.  Remarkably,  in  this  example  the  intuitive  idea  of 
displacing  parameter  Pi  away  from  the  stability  boundary  leads  to  the  loss  of  robustness  in 
C2  and  C3  because  the  tuning  adjustments  disregard  the  parametric  uncertainties  present. 

From  Theorem  2.2  it  follows  that  a controller  may  be  considered  optimally  robust  if  it 
is  capable  of  producing  the  smallest  possible  value  of  k n-  Hence,  the  problem  of  robust 
synthesis  reduces  to  the  following  optimization: 

c*  = argmin  &n(c)  (2.48) 

ce§° 

where  c*  is  the  optimal  controller,  J°  represents  the  set  of  all  controllers  which  lead  to 
stable  loops  with  the  nominal  plant,  and  &n(c)  is  the  robustness  parameter  (2.41)  for  a 
given  controller  c. 

The  constrained  optimization  problem  (2.48)  is  nonlinear  and  nondifferentiable.  In 
this  study  the  synthesis  problems  are  solved  by  exhaustive  numerical  search  over  the 
control  space.  This  approach  is  viable  only  because  of  the  low  dimensionality  of  the 
control  space.  For  higher  dimensions  this  approach  would  rapidly  become  computationally 
intractable  due  to  the  combinatorial  explosion  in  the  number  of  grid  points.  The  exhaustive 
search  is  carried  out  by  discretizing  the  space  of  control  parameters  with  a grid  of  size 
A(3o=A(3 1=0.5,  and  then  calculating  the  values  of  k^(c)  at  every  grid  point.  Grid  points 
corresponding  to  controllers  that  are  not  nominally  stable  are  discarded,  thereby  gaining 
execution  speed. 

The  last  row  of  Table  2.1  shows  the  numerically-determined  optimal  values  for  the 
control  and  robustness  parameters  for  design  C*.  The  optimal  design  is  robust  since  it 
satisfies  the  constraint  < 1,  and  it  achieves  the  smallest  value  of  &n  (0.501).  The 
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optimal  design  stabilizes  the  nominal  and  actual  plants,  as  is  confirmed  from  the  step- 
responses  shown  in  Figures  2.6a  and  2.6b. 


Table  2.1.  Control  and  robustness  parameters  for  conventional  controllers 
Cl,  C2,  and  C3,  and  optimally  robust  controller  C*. 


Controller 

Controller 

parameters 

Stability  margin 

Design 

P. 

Po 

% 

Cl 

101.0 

-98.0 

0.842 

C2 

115.0 

-98.0 

1.260 

C3 

127.0 

-98.0 

2.041 

C* 

31.0 

-30.5 

0.501 

Figure  2.7  shows  several  contour  plots  of  the  functional  &n(c)  as  a function  of  the 
controller  parameters  (3o  and  Pi.  The  contour  &n(c)  = 1 defines  the  boundary  of  robust 
stability.  The  optimal  design  is  remarkably  close  to  the  nominal  stability  boundary,  making 
it  very  sensitive  to  perturbations  in  the  control  parameters.  For  this  reason,  it  may  not  be 
desirable  to  design  a controller  to  optimize  the  stability  margin  alone. 


2.4.3  Discussion 

The  performance  in  servo-response  tests  is  documented  in  Figures  2.6a  and  2.6b  for 
the  various  designs  studied.  Only  the  nominal  and  actual  plants  are  considered.  The  figure 
shows  that  controller  C*  produces  low  overshoot  levels  and  minor  transient  oscillations. 
However,  caution  must  be  exercised  when  attempting  to  generalize  performance 
observations  from  this  example  because  the  optimally  robust  controller  has  been  designed 
ignoring  performance  constraints. 
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It  is  important  to  remark  that  functional  &n(c)  is  not  convex  with  respect  to  the  control 
vector  c.  Profiles  calculated  using  very  fine  grids  also  reveal  that  the  functionals  are  non 
smooth.  As  a consequence,  gradient-based  optimization  techniques  are  not  suitable 
substitutes  for  the  exhaustive  search  method. 


Figure  2.7.  Contour  plot  on  the  (Po,  Pi)  plane  of  the  functional  k^{c). 

2.5  Conclusions 

This  chapter  presents  a necessary  and  sufficient  condition  for  the  robust  stability  of  a 
continuous  time  system  subject  to  ellipsoidal  uncertainty.  The  condition  is  similar  to  that 
given  in  Guzzella  et  al.  (1991),  for  discrete  time  systems.  The  test  requires  a frequency 
search  to  be  performed.  In  order  to  avoid  the  problems  associated  with  such  a search,  a 
new  analysis  method  that  does  not  require  a frequency  search  is  presented  in  Chapter  3. 


CHAPTER  3 

STABILITY  MARGIN  CALCULATION  FOR  ELLIPSOIDAL  SYSTEMS 

3.1  Introduction 

The  goal  of  this  chapter  is  the  development  of  a robust  stability  test  for  ellipsoidally 
uncertain  systems  that  does  not  require  a frequency  search.  It  is  shown  that  the  necessary 
and  sufficient  condition  of  Theorem  2.1  can  be  transformed  into  a condition  of  the  form 

r«L  < i 

along  with  an  auxiliary  condition  for  co  = 0.  The  function  F(s ) is  a stable,  minimum- 
phase,  strictly  proper  transfer  function.  The  infinity  norm  condition  requires  only  the 
checking  of  the  eigenvalues  of  a certain  Hamiltonian  matrix,  and  thus  avoids  a frequency 
sweep.  Furthermore,  computing  the  actual  value  of  the  norm  can  be  done  using  a bisection 
search,  thereby  guaranteeing  a prespecified  level  of  accuracy.  Section  3.2  describes  the 
construction  of  the  transfer  function  F(s).  The  new  method  of  robust  stability  analysis 
using  this  transfer  function  is  developed  Section  3.3.  The  formulation  of  the  appropriate 
transfer  function  F(z)  for  discrete  systems  is  given  in  Section  3.4.  An  example  is  given  in 
Section  3.5  to  illustrate  the  proposed  method. 

3.2  Construction  of  E(s) 

To  test  the  necessary  and  sufficient  condition  for  robust  stability  established  in 
Theorem  2.1,  the  function  &N(oo)  must  be  evaluated  over  the  frequency  range  [0,°°).  This 
frequency  sweep  is  computationally  expensive  and  is  inherently  imprecise  since  a finite 
number  of  frequency  points  must  be  used.  An  alternative  is  to  construct  a stable  transfer 
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function  whose  frequency  response  is  equal  to  the  function  kN  (to)  for  positive  frequencies. 
Then  the  maximum  of  kN  (to)  for  positive  frequencies  can  be  found  as  the  infinity  norm  of 
this  transfer  function.  Actually,  as  a first  step  in  this  construction,  the  square  of  the 
function  &N(©)  will  be  fitted  to  make  the  mathematical  developments  more  tractable.  From 


the  resulting  form,  the  appropriate  F(s ) can  be  easily  constructed. 
The  function  &N(w)  is  defined  in  Chapter  2 as 

1 


*n(®)  = 


XT(0))Qco1X(03) 


a/  tfk.,1 


So 


© > 0 
© = 0 


where 


x(©)  = 


'Re  {G°0'©)}‘ 

'X/j(  ©)' 

Im  {G0  (;©)} 

X/(©) 

e W 


(3.1) 


(3.2) 


and 


A(©) 

B(©) 

wJ(©)Qgw^(©) 

wJ(©)QgW/(©) 

B(w) 

D(©)_ 

w]’(w)Qgw^(©) 

w}  (©)Qgw/(©) 

(3.3) 


The  vectors  w/?(©)  and  wX©)  are  defined  in  (2.18). 

Define  /z(©)  to  be  the  square  of  the  function  that  defines  &N(©)  for  positive 
frequencies 


M©)  = ~r ~\ 

xT«o)Q(01t(©) 


(3.4) 


Even  though  Q'1  is  not  defined  at  © = 0,  a limiting  value  of  /i(©)  as  © ->  0 is  well 
defined,  but  this  limiting  value  is  not  equal  to  &N(0)  in  general.  Substituting  equations 
(3.2)  and  (3.3)  into  (3.4)  yields 


h((0)  = 


A(©)D(©)-B2(©) 

A(©)t2  (©)  + D(©)x^(©)  - 2B(©)x/?(©)t/(©) 


(3.5) 


The  objective  is  to  find  a stable  transfer  function  F(s ) that  satisfies 
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|F(7-co)|2  = h(  co)  (3.6) 

for  non-zero  frequencies.  To  facilitate  construction  of  such  an  F(s),  Lemma  3.1  states 
useful  properties  of  the  function  h( co). 

Lemma  3.1  The  function  h(c o)  is  a finite,  non-negative,  even,  rational  function  of 
(0,  with  numerator  degree  4k-8  and  denominator  degree  4k-4. 


Proof:  The  proof  of  Lemma  3. 1 is  given  in  the  Appendix.  V 

Lemma  3.1  is  used  in  the  proof  of  Theorem  3.1,  which  describes  a factorization  of  h{ to) 
used  to  compute  F(s). 


Theorem  3.1  The  function  //({ o)  has  a factorization  of  the  form 

*»)  - 


mv 


(3.7) 


where  the  polynomials 

mx 

X(s)  = ^xrsr  and  Y(s)  = ^ yr/ 

r = 0 r = 0 

have  no  zeros  in  the  closed  right  half  plane  (RHP),  have  real  coefficients,  and  mx  = 2.k-4 
and  mY  = 2k-2. 


Proof:  The  proof  is  based  on  a result  from  spectral  factorization,  given  by  Rozanov 
(1967),  which  states  that  every  non-negative,  even,  rational  function  of  a scalar  variable  co 
has  a factorization  of  the  form  (3.7)  where  X(s)  and  Y(s)  have  no  zeros  in  the  open  RHP 
and  have  real  coefficients.  Lemma  3.1  shows  that  h{ co)  satisfies  the  conditions  of 
Rozanov,  and  thus  has  such  a factorization.  The  degree  of  the  numerator  of  h{ co)  is  4&-8 
and  this  is  twice  the  degree  of  X(s).  Similarly,  the  denominator  degree  of  h( co)  is  4k-4, 
and  this  is  twice  the  degree  of  Y(s).  Furthermore,  since  h( co)  is  finite  and  non-zero  for  all 
co,  it  follows  that  X(s)  and  Y(s)  have  no  zeros  in  the  closed  RHP.  V 


If  the  function  F(s ) is  constructed  as 
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then 


F(s) 


X(s) 
Y (s) 


(3.8) 


\F(jaf 


X(yto)X(-y‘(o) 

X(y'co) 

Y(yto)Y(-yto) 

YO'co) 

= h(c o) 


(3.9) 


and  (3.6)  is  satisfied.  It  remains  to  calculate  the  coefficients  of  the  polynomials  X(s)  and 


Y(s)  so  that  (3.7)  is  satisfied.  The  first  step  is  to  express  h( to)  as  a ratio  of  polynomials  in 
yco  so  that  the  coefficients  of  X(/( o)  and  Y(/G)),  and  thus  of  X(s ) and  Y(s),  can  be  found. 


Theorem  3.2  The  function  /i(oo)  defined  in  (3.5)  is  equivalent  to 


h{  co)  = — 


A(;co)DO'co)-B20-co) 


A(y(0)T/(yco)  + D(;0)x^(y(O)  - 2B(;o))x^(;co)x/(y'co) 


(3.10) 


where  the  coefficients  of  the  polynomials  A(y'(0),  B(y'co),  D(yc o),  and  X/0'(fl) 

depend  only  on  the  nominal  coefficient  vector  g°  and  the  matrix  Qg.  The  polynomials 
XflO'G))  and  X/Oco)  are  given  by 


mi 


* rU °))  - Xs^Ca0)2 

(=0 


with  mR  = < 


-k  k even 

, 2 (3.11a) 

±(k- 1)  k odd 


and 


W®)  = Xg?2£-l)0'co) 
( =1 


with  in  | = 


-k  A:  eve  n 

. 2 (3.11b) 

.^■(fc  + l)  A:  odd 


where  g,°  are  the  elements  of  the  vector  g°  defined  in  (2.15),  except  for  g®  :=  1.  The 
polynomials  A(;co),  B(y'co),  and  D (joo)  are  given  by 


and 


A(;(0)  = Xa2^0'“) 
(=0 


2( 


with  n a = 


\k-  2 A:  even 
I k - 1 k odd 


(3.12a) 


and 


BO'co)  = Xb(2M)0'“) 

e=\ 


.v.^(2^-i) 


D(yto)  = £d2  eU<°) 
e=\ 


2( 


with  nB  = k - 1 


with  nD  = 


[ k - 1 k even 
\k- 2 A:  odd 


(3.12b) 


(3.12c) 
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Furthermore,  the  coefficients  of  the  polynomials  A(y'co),  B(y'co),  and  6(703)  are  given  by 


a2(  - • 


^(2M)  “ < 


X 9(k-2h),(k-2£+2h) 

h=0 

mA 

X Cl(k-2h),(k-2f+2h) 

h=OmA) 


X Cl(k-2h),(k-2^+2h+l) 

h=0 

m2 

X Cl(k-2h),(k-2^+2h+l) 

h=(£-m]) 


0 < i < mA 
(mA  + l)<^<nA 

1 < ^ < m, 

(irq  + 1)  < l < nB 


(3.13a) 


(3.13b) 


d2^  - 


1 

X (l(k-2h+l),(k-2^+2h-l) 
h=l 

mD 

X tl(k-2h+l),(k-2^+2h-l) 


1 < t < mD 
(mD  + l)<^<nD 


(3.13c) 


h=(^-mD+l) 

respectively,  where  qjj  represents  the  (i,j)  element  of  the  matrix  Qg,  and  the  indices  mA, 
mo,  mi  and  m2  are  given  by 


and 


mA  = 2nA 


mD  = 2<nD  + 1) 


m>  = 


(2 


(k- 1) 


m2  = 


±(k-  2) 


k even 
k odd 

k even 
k odd 


Proof:  The  proof  of  Theorem  3.2  is  given  in  the  Appendix. 

Theorem  3.2  gives  an  expression  that  can  be  used  to  compute  the  function  h{ (0)  as  a 
rational  function  of 7(1).  However,  the  coefficients  of  the  polynomials  X(s)  and  Y(s)  are  the 
desired  quantities.  It  is  not  possible  to  extract  a simple  expression  for  the  coefficients  of 
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X(s)  and  Y(s)  in  terms  of  the  variable  c,  however.  Therefore,  the  coefficients  of  these 
polynomials  must  be  found  numerically. 

The  problem  of  computing  the  coefficients  of  the  stable  polynomials  X(s)  and  Y(s) 
given  the  coefficients  of  the  polynomials  |X(y(D)|  = X(j(0)X(-j(0)  and 
|Y(yO))|  = Y(yO))Y(-yO))  is  a spectral  factorization  problem.  Efficient  algorithms  for 
solving  this  problem,  for  discrete  or  continuous  systems,  are  given  in  Jezek  and  Kucera 
(1985)  and  Kucera  (1979).  The  algorithms  presented  in  Jezek  and  Kucera  (1985)  are 
iterative,  but  are  quadratically  convergent,  and  can  be  used  to  compute  the  desired 
polynomial  coefficients  to  a user-specified  tolerance. 


3.3  Robust  Stability  Testing  Using  F(s ) 


It  follows  from  Section  3.2  that  F(s)  defined  by  equation  (3.8)  is  strictly  proper,  stable 
and  minimum  phase.  Now,  since  h( CD)  is  taken  to  be  the  square  of  the  function  that  defines 
x(co),  it  follows  that 

XU©) 


\F(M  = 


YU©) 

for  nonzero  frequencies,  implying  that 


= X(;0))X(  •/(°)  = VW)  = *(©)  (3-14) 

i YU©)Y(-y©) 


sup  *(©)  = \\F(s)\\o 
co>0 


(3.15) 


Equations  (3.15)  and  (3.1)  for  to  = 0 imply  that 


&N  = max  &N(co)  = max 
co 


linolL . 


-Jtfk.k 

bo  I 


(3.16) 


so  that  testing  the  condition  of  Theorem  2. 1 is  equivalent  to  testing  the  two  conditions 

\\F(s)\l  < 1 (3.17a) 


So 


< 1 


(3.17b) 
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The  function  F(s ) is  a strictly  proper,  stable  transfer  function,  so  the  algorithm 
proposed  by  Boyd  et  al.  (1989)  can  be  used  to  compute  ||ir(s)||oo.  This  algorithm  requires  a 
state-space  realization  of  F(s ) as  input.  The  basic  idea  behind  the  algorithm  is  that  if  a 
strictly  proper,  stable,  real-rational  transfer  function  H(s)  has  the  realization 

H(s)  = C(sl-\y'B  (3.18) 

then 


ItfWlL  S Y 


if  and  only  if  the  matrix 


My  := 


y 'BB1 


_-y  C C -A' 

has  at  least  one  purely  imaginary  eigenvalue.  This  implies  ||//(.s)||oo  < y if  and  only  if  My 
has  no  purely  imaginary  eigenvalues.  The  function  F(s)  will  have  a state-space  realization 
of  the  form  given  in  equation  (3.18)  since  it  is  strictly  proper.  F(s)  is  stable  also,  so  the 
algorithm  can  be  used  to  test  the  condition  (3.17a)  by  checking  if  the  matrix 


M := 


A BB1 
-CTC  -A1 


(3.19) 


has  any  purely  imaginary  eigenvalues.  The  actual  value  of  ||F(s)|L  IS  f°und  through  a 
bisection  search  on  y.  Upper  and  lower  bounds  for  ||F(s)|L  in  terms  of  the  realization 
{A,B,C}  are  given  in  Boyd  et  al.  (1989).  A method  for  testing  whether  or  not  My  has 
any  imaginary  eigenvalues  without  actually  computing  the  eigenvalues  is  also  given.  Such 
a method  could  be  used  to  avoid  any  numerical  problems  computing  eigenvalues.  The 
robust  stability  test  that  has  been  developed  is  summarized  below. 

Robust  Stability  Testing  Algorithm 

1)  Obtain  the  nominal  plant,  the  uncertainty  matrix  Q^,  and  a nominally  stabilizing 
controller  C(s). 


2)  Construct  the  vector  g°  and  the  matrix  Qg. 
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3)  Compute  the  coefficients  of  the  polynomials  A(yco),  B(y(o),  D(yco),  x^(y'co)  and 
T/O'co)  given  in  (3.11  )-(3. 1 3). 

4)  Construct  the  function  h{ to)  given  in  equation  (3. 10). 

5)  Perform  spectral  factorizations  to  find  Xfs)  and  Y(s)  with  roots  in  the  left  half  of  the 
complex  plane. 

6)  Construct  F(s)  and  obtain  {A,B,C}  such  that  F(s)  = CfsI-Af'B. 

7)  Calculate  iFfsJL  using  the  algorithm  of  Boyd  et  al.  ( 1989). 

8)  Robust  Stability  Test:  Test  the  conditions  given  in  (3.17a)  and  (3.17b). 

The  advantage  of  this  method  for  testing  robust  stability  is  that  it  avoids  having  to 
construct  a frequency  plot  or  perform  a grid  search.  Testing  condition  (3.17a)  requires 
only  checking  the  eigenvalues  of  one  Hamiltonian  matrix  constructed  from  the  state-space 
matrices  of  F(s ).  Calculating  the  coefficients  of  the  polynomials  used  to  find  F(s)  is 
straightforward.  Furthermore,  computing  the  actual  value  of  the  norm  in  (3.17a)  can  be 
done  using  a bisection  search  based  on  the  algorithm  mentioned  above.  Thus,  the  stability 
margin  of  the  system  can  be  obtained  to  any  desired  degree  of  accuracy. 


The  necessary  and  sufficient  condition  for  robust  stability  for  the  discrete  time  case  is 


where  pc(co)  is  given  by  equation  (2.37).  This  condition  is  equivalent  to  the  three 
conditions 


3.4  Discrete  Time  Stability  Analysis 


pc(to)  < G°(eju>)  V to  e [0,tc] 


(3.20) 


< 1 


(3.21a) 


< 1 


M0)| 


(3.21b) 
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M*)  I 


< 1 


(3.21c) 


As  in  the  continuous  time  case,  the  endpoint  conditions  (3.21b)  and  (3.21c)  will  be  left  to 
test  separately.  Defining  h( CO)  as 


h(  to)  = — y 

^ (o))Q(01x(to) 

then  it  follows  from  the  definitions  in  Chapter  2 that 


(3.22) 


h{  (0)  = 


A(co)D(co)-B  (co) 


A(co)i/  (co)  + D(co)x|(co)  - 2B(co)x/?(co)x/(co) 


(3.23) 


The  goal  is  to  construct  a stable  transfer  function  F(z)  that  satisfies 

i2 


\F(ej(a)\  = /i(co)  (3.24) 

Invoking  the  discrete  time  version  of  Theorem  3.1  given  in  Rozanov  (1967),  it  follows  that 
/i(co)  has  a factorization  of  the  form 


K co)  = 


_ [X(^M)| 

|Y(<>)| 


(3.25) 


where  the  polynomials 


mx 


IK'S 


X(z)  = ^ and  Y(z)  = ^yr  zr 

r = 0 r = 0 

have  all  their  zeros  in  the  open  unit  circle.  The  problem  is  treated  in  a manner  similar  to  that 
given  in  Section  3.2,  except  now  the  function  h{ co)  is  rewritten  as  a function  of  c^“. 


Theorem  3.3  The  function  h{ co)  defined  in  (3.22)  is  equivalent  to 

A(ejti>)D(ej<0)  - B2(ej(i>) 


h(  co)  = 


(3.26) 


A(eJm  )x]  (ej(0 ) + D(ej(a  )x2R  (ej(0 ) - 2B  (ejl “ )x  R (ejm  )x , (eja ) 
where  the  coefficients  of  the  polynomials  A(e ■yco),  B(e7C0),  D(e-/t0),  xR(e^)  and  xt{e^) 
depend  only  on  the  nominal  coefficient  vector  g°  and  the  matrix  Qg.  The  polynomials 
xR(e^m)  and  x/(e'/t0)  are  given  by 
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M*yco) = + *~y£(0)  0.27a) 

£=0 

and 

T/^05)  = £g?(^“  - e~Jea)  (3.27b) 

^=o 

where  g°  are  the  elements  of  the  vector  g°  defined  in  (2.13),  except  for  g®  :=  1.  The 
polynomials  A(e;(0),  B(e^),  and  D(eya))  are  given  by 


and 


and 


k- 1 2k-2 

= 4qu  + 2Xq„  + !>,(<>“  + e'*") 
£=1  r=l 

2k-2 

B(eja)  = Xbr(^'r(0  - e~jm) 

r=l 

k- 1 2*-2 

D(^)  = -2£q„  + Xdr(^“  + «->“) 

*=1  r=l 


(3.28a) 


(3.28b) 


(3.28c) 


The  coefficients  ar,  br,  and  dr  are  determined  by  the  magnitude  of  r.  Define  the  quantities 


a 


k 

X292*-r-« 

^=k-(m  f -1 ) 


V 


k- r 

= X2cU  r+( 

e=i 


(3.29) 


where  mj  = (r/2)  if  r is  even  and  mi  = (r- 1 )/2  if  r is  odd.  Now,  if  1 < r < k-l,  then  the 
coefficients  are  given  as 


H 

| 9k-m2,k-m2 

[ G + V 

+ <J  + V 

r even 
r odd 

(3.30a) 

b,  = 

J 9k-/n2,k-m2 

\<r 

+ a 

r even 
r odd 

(3.30b) 

9k-m2,k-m2 

o - o 

+ 0-0 

r even 
r odd 

(3.30c) 

Finally,  if  k < r < 2k-2,  then 
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ar  = br  = dr  = 


L 

Qrri2,ni2  + X2ck2k-r- 


1^2  ,A722 
m2 


X2ck2k-r-f 

U=1 


r even 


r odd 


(3.31) 


where  m2  = k - (r/2)  if  r is  even  and  m2  = k - (r+l)/2  if  r is  odd. 
Proof:  The  proof  of  Theorem  3.3  is  given  in  the  Appendix. 

When  h{ co)  is  computed  according  to  equation  (3.26),  the  result  is 


h{  co)  = 


4*-4 

X"< 

e=o 

ej(ae  + e-W' 

e-j(4k-4)(0< 

r4k-4 
J= 0 

' ej((+4k-4)m  + g-;«-4A:+4)(oj| 

4k-2 

x^ 

£=0 

'ej(oe  + e~W\ 

e-j(4k-2)(o< 

4k-2 

Id, 

J= 0 

' ej((+4k-2)(0  + e-j((-4k+2)w j , 

where  n(  and  d^  are  the  coefficients  that  result  from  equations  (3.27)-(3.31).  This  is 
equivalent  to 
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for  appropriately  calculated  nq  and  d(.  As  a result  of  the  fact  that  many  of  the  polynomial 
coefficients  given  by  (3.28)-(3.31)  are  identical  up  to  sign,  many  of  the  coefficients  in 
(3.32)  are  exactly  zero.  The  number  of  terms  lost  is  the  same  for  both  the  numerator  and 
denominator,  and  is  equal  to  4(/c-l).  Removing  these  terms  gives  a reduced  form  for  h( co) 

(3.33) 
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The  polynomials  X(e/C0)  and  Y(e/C0)  are  found  by  performing  a spectral  factorization  using 
the  results  of  Jezek  and  Kucera  (1985)  in  a manner  similar  to  the  continuous  time  case. 
The  final  polynomials  X(z)  and  Y (z)  have  degree  2k- 1,  however,  the  first  coefficient  of 
XU)  is  identically  zero. 
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The  transfer  function  that  results  in  the  discrete  case  is  of  comparable  order  to  that 
obtained  in  the  continuous  case,  and  this  order  is  roughly  twice  that  of  the  nominal 
characteristic  polynomial.  In  both  cases,  the  transfer  function  is  found  by  constructing  a 
stable  transfer  function  from  the  vector  of  coefficients  of  the  nominal  characteristic 
polynomial  and  the  elements  of  the  matrix  describing  the  uncertainty  ellipsoid. 


3.5  Example 

In  this  section,  an  example  is  presented  to  illustrate  the  ideas  introduced  in  this  chapter. 
Consider  the  following  nominal  polynomial 

G°(s)  = s4  + 5.8600s3  + 9.3954s2  + 6.0126s  + 5.3237  (3.34) 

represented  by  g°,  and  the  associated  uncertainty  matrix  Qg: 


"5.8600' 

"2.0425 

2.3648 

1.7252 

1.2603" 

g°= 

9.3954 

2.3648 

4.9454 

3.5362 

1.4123 

6.0126 

5.3237 

Qg  = 

1.7252 

1.2603 

3.5362 

1.4123 

4.6202 

2.2330 

2.2330 

3.6206 

The  polynomial  (3.34)  has  roots  at  s = {-3.60,  -2.00,  -0. 13±0.85y}.  Using  the 
equations  from  Theorem  3.2  and  the  algorithm  of  Jezek  and  Kucera  (1985),  the  spectral 
factors  X(s)  and  Y(s)  for  this  system  are  found  to  be 

X(s)  = 2.12s4  + 6.37s3  + 9.58s2  + 7.84s + 3.43  (3.35) 

and 

Y(s)  = 1. 43s6  + 6.92s5  + 2 1. 7s4  + 34.2s3  + 34.9s2  + 23.3s + 10.9  (3.36) 

The  degree  of  the  nominal  polynomial  in  equation  (3.34)  is  k = 4;  it  is  easily  verified  that 
the  degrees  of  X(s)  and  Y(s)  match  the  specifications  of  Theorem  3.1.  The  roots  of  X(s) 
are  s = {-1 ,02±0.60/,  -0.48±0.96y},  and  those  of  Y(s)  are  s = {-1 .30±2.22y, 
-0.99±0.63/,  — 0. 13±0.90/} , showing  that  both  polynomials  are  Hurwitz. 
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Figures  3.1a  and  3.1b  show  plots  of  the  frequency-dependent  Nyquist  robust  stability 
margin  vs.  frequency  for  the  frequency  range  CD  e [0,  10].  It  is  obvious  that  k^((0)  is 
discontinuous  at  CD  = 0 for  this  example,  and  k^(0)  = 0.357420.  Using  equation  (3.16), 
the  value  of  &n  is  found  to  be  k n = 0.971650,  rounded  to  six  digits,  and  since  the 
condition  (3.17)  is  satisfied,  the  polynomial  is  robustly  stable.  The  plot  in  Figure  3.1a  is 
constructed  by  evaluating  (3.1)  over  the  frequency  range  CD  e [0,  10]  using  Acd  = 0.01. 
The  plot  of  the  frequency  response  magnitude  of  X(s)/Y(s),  evaluated  at  the  same 
frequency  points,  is  shown  in  Figure  3.1b.  Figures  3.1c  and  3. Id  show  a comparison  of 
the  results  of  the  two  methods.  The  maximum  absolute  difference  between  the  two  curves 
is  less  than  4xl0'6,  while  the  maximum  percentage  difference  is  slightly  greater  than 
0.0004%. 

The  results  of  the  frequency  sweep  method  are  shown  in  Table  3.1  for  different  values 
of  Acd.  The  computed  value  of  k n (rounded  to  six  digits)  as  well  as  the  number  of  floating 
point  operations  needed  by  a simple  MATLAB  routine  are  shown.  It  is  interesting  to  note 
that  the  result  for  Acd  = 0.3  is  quite  good,  and  in  fact,  much  better  than  the  result  for 
Acd  = 0.2.  This  is  due  to  the  fact  that  the  maximum  of  &n(cd)  occurs  almost  exactly  at 
CD  = 0.9.  This  point  will  be  evaluated  in  the  sweep  using  Acd  = 0.3,  yielding  a much 
better  estimate  of  k n than  for  the  sweep  using  Acd  = 0.2. 


Table  3.1.  Floating  point  computations  needed  for  a frequency  search. 


Resolution  Acd  for  the  frequency  sweep 

0.3 

0.2 

0.1 

0.01 

0.971648 

0.787774 

0.971648 

0.971648 

flops 

14,425 

21,858 

43,718 

437,206 

abs.  error  % error 
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Figure  3.1.  Plot  of  the  Nyquist  robust  stability  margin  &n(0))  vs.  frequency  to. 

(a)  Frequency  sweep  method,  (b)  Proposed  method,  (c)  Percentage 
error  between  the  two  methods,  (d)  Absolute  error  between  the  two 
methods. 
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Tables  3.2  and  3.3  show  the  results  for  the  method  proposed  in  this  paper.  In  Table  3.2, 
the  number  of  floating  point  computations  needed  to  find  both  spectral  factors  is  listed. 
These  numbers  represent  the  computations  of  a simple  MATLAB  program  implementing 
the  algorithm  given  in  Jezek  and  Kucera  (1985).  The  tolerance  is  taken  to  be  the  sum  of 
squared  errors  between  the  generated  coefficients  of  X(yco)X(-yco)  or  Y(yco)Y(-_/co)  and 
the  actual  values  of  the  coefficients.  This  tolerance  must  be  satisfied  by  both  polynomials 
X(s)  and  Y(s).  The  number  of  iterations  of  the  algorithm  in  Jezek  and  Kucera  (1985) 
required  to  attain  this  tolerance  is  also  given  in  Table  3.2.  For  the  computations  presented 
in  Table  3.3,  the  value  IO-4  was  used  for  the  tolerance. 

Table  3.2.  Floating  point  computations  and  iterations  needed  to  compute  X(s)  and  Y(s). 


Tolerance  for  computing  spectral  factors 

10-2 

io-4 

io-6 

io-9 

iterations  (X/Y) 

4/6 

4/6 

5/7 

5/7 

flops 

3682 

3682 

4337 

4337 

Table  3.3  shows  the  number  of  floating  point  computations  needed  to  find  IIF(.s)IL 
given  X(s)  and  Y(s).  The  tolerance  values  listed  are  the  relative  tolerances  required  for 
termination  of  the  algorithm  used  to  compute  IIF(5)IL,  i.e.,  the  estimated  value  of  IIF(.y)IL  is 
guaranteed  to  be  within  ±£llF(s)IL  of  the  true  value.  It  is  very  interesting  to  note  that  the 
number  of  computations  required  to  find  the  spectral  factors  X(s)  and  Y(s)  is  small  relative 
to  the  computation  required  to  find  IIF(5)IL.  This  is  desirable,  as  the  termination  criterion 
for  the  spectral  factorization  algorithm  can  be  stringent  without  seriously  affecting  the 
number  of  overall  calculations  needed  to  compute  IIF(.s)IL.  Thus,  the  computational  burden 
of  the  proposed  method  is  primarily  dictated  by  the  desired  numerical  accuracy  of  the 
Nyquist  robust  stability  margin. 
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Table  3.3.  Floating  point  computations  needed  for  computing  IIF(s)IL. 


Relative  tolerance  £ for  computing  IIF(s)IL 

0.1 

0.01 

0.001 

0.00001 

WF(s)\L 

0.917423 

0.973195 

0.971801 

0.971659 

flops 

181,669 

268,629 

358,405 

573,145 

For  a tolerance  as  low  as  e = 0.001  used  in  computing  the  infinity  norm,  the 
computational  burden  for  the  proposed  algorithm  (358,405  + 3,682  = 362,087  flops)  is 
lower  than  that  for  the  frequency  sweep  using  Aco  = 0.01  (437,206  flops).  If  the  relative 
tolerance  level  is  increased  to  e = 0.00001,  then  the  computational  burden  (573,145  + 
3,682  = 576,827  flops)  surpasses  that  of  the  frequency  sweep  using  Aco  = 0.01. 
However,  the  proposed  algorithm  affords  a guaranteed  level  of  accuracy  in  the  final  result 
that  the  frequency  sweep  cannot. 

3.6  Conclusions 

In  this  chapter,  a new  robust  stability  test  for  ellipsoidal  systems  is  presented.  This 
test  is  based  on  construction  of  a strictly-proper,  stable  transfer  function  whose  infinity 
norm  is  equivalent  to  the  stability  margin  of  the  system.  The  proposed  test  avoids 
performing  a frequency  search,  and  can  calculate  the  stability  margin  using  a bisection 
search. 


CHAPTER  4 

ROBUST  PREDICTIVE  CONTROL  DESIGN  FOR  ELLIPSOIDALLY 

UNCERTAIN  SYSTEMS 


4.1  Introduction 

The  issue  of  designing  controllers  to  robustly  stabilize  a system  subject  to  real 
parametric  uncertainty  has  received  much  attention  in  recent  years.  Since  many  physical 
uncertainties  are  most  directly  and  effectively  modeled  by  real  parametric  uncertainty,  it  is 
of  great  practical  interest  to  include  real  uncertainty  in  the  design  of  robust  controllers. 
Unfortunately,  the  inclusion  of  real  uncertainties  greatly  complicates  both  the  robust  control 
analysis  and  synthesis  procedures.  However,  there  are  cases  of  interest  where  the 
uncertainty  structure  is  such  that  the  analysis  and  synthesis  problems  simplify 
considerably.  This  chapter  considers  such  a case  and  establishes  a procedure  to  design  a 
robust  predictive  controller  for  a system  subject  to  ellipsoidal  parametric  uncertainty. 

The  standard  method  for  analyzing  the  effect  of  uncertainty  on  a system  is  to  construct 
a fictitious  feedback  loop  where  all  the  uncertainties  are  pulled  out  in  a block-diagonal  A 
matrix  and  all  the  known  components,  including  the  nominal  plant,  the  controller,  and  any 
weights,  are  lumped  into  an  M matrix.  The  robust  control  analysis  problem  is,  given  a 
fixed  M and  a set  structure  for  A,  find  the  smallest  A that  destabilizes  the  system.  The 
structured  singular  value  (Doyle  1982),  p,  is  a stability  margin  defined  by  the  size  of  the 
smallest  destabilizing  A as 

p(M)  :=  {min  a(A)  I det(I  - AM)  = 0}-1 

where  p(M)  = 0 if  no  allowable  A destabilizes  the  system.  The  robust  control  synthesis 
problem  is,  given  a description  of  all  possible  As,  design  a controller  that  stabilizes  the 
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closed  loop  for  any  allowable  A.  It  is  the  structure  of  the  A matrix  that  determines  the 
difficulty  associated  with  solving  the  analysis  and  synthesis  problems. 

The  simplest  uncertainty  structure  is  where  the  matrix  A is  a full  complex  block,  and 
only  a norm  bound  on  the  whole  matrix  is  known.  This  case  is  referred  to  as  unstructured 
uncertainty,  since  no  internal  structure  for  the  A matrix  is  known  (or  assumed).  In  this 
case,  both  the  analysis  and  synthesis  problems  simplify  considerably.  The  structured 
singular  value  for  a full  complex  block  A is  simply  p(M)  = a(M).  Furthermore,  a 
controller  is  robustly  stable  if  it  generates  an  M matrix  that  satisfies  HMf^  < flAf”1 . Thus, 
Hoo  design  methods  (Doyle  et  al.  1989),  (Glover  and  Doyle  1988)  can  be  used  to  find  the 
optimal  controller,  or  any  sub-optimal  controller  that  stabilizes  the  closed  loop. 

In  many  cases  of  practical  interest,  the  block-diagonal  A matrix  contains  real  and 
complex  (possibly  repeated)  scalar  blocks  as  well  as  norm  bounded  complex  uncertainty 
blocks.  If  A contains  real  blocks  as  well  as  complex  blocks,  p(M)  is  referred  to  as  mixed 
p.  In  general,  the  calculation  of  mixed  p is  NP  hard  (Braatz  et  al.  1994),  (Braatz  1996), 
(Poljak  and  Rohn  1993).  For  both  the  purely  complex  and  mixed  cases,  convex  upper 
bounds  for  p exist  (Doyle  1982),  (Fan  et  al.  1991),  and  these  are  the  basis  for  controller 
synthesis  methods.  The  complex  p-synthesis  control  design  procedure  attempts  to  find  a 
controller  minimizing  p for  the  closed  loop  by  iterating  between  computing  an  optimal 
scaling  matrix  D appearing  in  the  upper  bound  and  computing  an  Hoo  controller  K for  this 
D,  hence  the  name  D-K  iteration  (Doyle  1985).  The  procedure  is  not  jointly  convex  in  D 
and  K so  there  is  no  guarantee  that  the  resulting  controller  is  globally  optimal,  but  it  will  be 
a locally  optimal  choice.  The  process  of  controller  design  for  mixed  p is  similar  to  that  for 
complex  p,  but  now  another  scaling  matrix  is  involved  in  the  optimization,  leading  to  the 
so-called  D,G-K  iteration  method  (Young  1996).  This  optimization  is  not  jointly  convex  in 
D,  G and  K,  so  it  may  not  yield  the  globally  optimal  controller. 

In  the  structured  singular  value  framework,  the  individual  elements  of  the  A matrix  are 
assumed  to  vary  independently.  Often  it  is  more  realistic  to  allow  interactions  between  the 
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uncertainties.  One  practical  example  is  the  ellipsoidal  uncertainty  description  where  the 
allowable  coefficients  of  the  system  transfer  function  are  restricted  to  lie  inside  an  ellipsoid 
in  the  plant  parameter  vector  space.  This  corresponds  to  a A matrix  with  real  diagonal 
elements  that  is  constrained  by  a weighted  two-norm  bound.  Ellipsoidal  uncertainty 
descriptions  arise  naturally  in  parameter  identification  techniques  where  the  uncertainty 
ellipsoid  is  associated  with  the  parameter-error  covariance  matrix.  Ellipsoidal  uncertainty 
descriptions  are  commonly  encountered  in  chemical  engineering  applications  where  model 
parameters  are  found  by  fitting  experimental  data  using  linear  or  nonlinear  regression 
techniques  (Fogel  and  Huang  1982),  (Belfonte  and  Bona  1985),  (Belfonte  et  al.  1990). 
Methods  for  the  computation  of  stability  margins  for  ellipsoidal  systems  are  given  by 
several  authors  such  as  Guzzella  et  al.  (1991),  Biernacki  et  al.  (1987),  and  Tsypkin  and 
Polyak  (1991).  Biernacki  et  al.  (1987)  outline  an  iterative  robust  controller  design 
procedure  based  on  expanding  the  parameter  space  stability  hyperellipsoid  until  it  covers 
the  allowable  uncertainty  region,  but  this  procedure  is  not  a convex  optimization.  The 
generalized  structured  singular  value  (Chen  et  al  1994a,  Chen  et  al  1994b)  is  an  extension 
of  the  structured  singular  value  that  takes  into  account  interactions  of  the  system 
uncertainties.  The  complex  uncertainties,  real  uncertainties  and  complex  block  uncertainties 
are  grouped  together,  and  each  group  is  norm  bounded  in  an  appropriate  manner. 
Therefore,  this  is  a more  general  framework  than  the  mixed  |i  description.  Recently, 
Braatz  and  Crisalle  (1997)  recast  the  ellipsoidal  uncertainty  problem  in  terms  of  the 
generalized  structured  singular  value,  allowing  the  straightforward  inclusion  of  complex 
uncertainties  in  the  stability  margin  calculation.  However,  there  are  no  currently  available 
controller  design  methods  based  on  the  generalized  structured  singular  value. 

In  some  specialized  cases,  such  as  when  the  uncertain  parameters  appear  in  an  affine 
manner  in  the  numerator  and  denominator  of  a single  input,  single  output  system  transfer 
function,  the  M matrix  that  results  is  rank  one.  Then  the  structured  singular  value  is  equal 
to  its  convex  upper  bound  (Young  1994),  (Young  and  Doyle  1996),  and  can  be  computed 
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directly.  The  generalized  structured  singular  value  can  also  be  computed  in  closed  form 
(Chen  et  al.  1994b),  (Braatz  and  Crisalle  1997).  It  is  important  to  note  that  the  restriction 
of  the  M matrix  to  be  rank  one  does  not  depend  on  the  number  of  uncertain  parameters,  but 
rather  the  manner  in  which  they  enter  the  system  description.  For  the  case  of  a rank  one  M 
matrix,  Rantzer  and  Megretski  (1994)  derive  a convex  necessary  and  sufficient  condition 
for  robust  stability  that  is  used  as  the  basis  for  design  of  a fixed-order  Youla  parameter 
whose  coefficients  are  found  using  convex  or  quasi-convex  optimization.  The  uncertainty 
description  considered  in  this  chapter  yields  a rank  one  M matrix,  and  the  method  of 
Rantzer  and  Megretski  (1994)  is  used  to  compute  an  appropriate  Youla  parameter  (Youla  et 
al.  1976)  that  defines  a robust  predictive  controller. 

Predictive  control  is  a model-based  control  methodology  that  has  found  wide 
acceptance  in  industry.  This  is  because  predicitive  control  offers  good  performance,  is 
easy  to  understand  and  formulate,  and  can  accommodate  input/output  process  constraints. 
The  industrial  success  of  the  predictive  control  techniques  is  apparent  by  the  variety  of 
commercial  predictive  controllers  that  are  available  to  the  chemical  processing  industry 
through  specialized  vendors.  Seborg  (1994)  reports  that  in  oil  refineries  and  petrochemical 
plants  around  the  world,  there  are  hundreds,  perhaps  thousands,  of  predictive  controllers 
employed. 

It  is  possible  to  design  predictive  controllers  using  different  plant  representations, 
including  finite  impulse  response  (FIR)  models,  transfer  function  models  and  state  space 
models.  FIR  based  schemes  include  Dynamic  Matrix  Control  (DMC)  (Cutler  and  Ramaker 
1980),  Model  Algorithmic  Control  (Mehra  et  al.  1979),  and  the  quadratic  DMC 
formulation  of  Garcia  and  Morshedi  (1986).  These  methods  are  applicable  only  to  stable 
plants,  and  often  require  large  model  orders,  typically  retaining  30  to  40  impulse  response 
coefficients.  Transfer  function  models  are  applicable  to  both  stable  and  unstable  plants, 
and  lead  to  lower  order  representations.  The  well-known  Generalized  Predictive  Control 
(GPC)  technique  (Clarke  et  al.  1987)  and  the  MUSMAR  approach  (Greco  et  al.  1984)  are 


53 


examples  of  transfer  function  based  predictive  control  formulations.  Kwon  and  Pearson 
(1977)  and  Muske  and  Rawlings  (1993)  present  state  space  formulations  for  predictive 
control.  This  large  body  of  literature  constitutes  a rich  source  of  knowledge  to  support  the 
design  and  analysis  of  predictive  controllers. 

There  is  currently  a large  amount  of  research  focusing  on  the  issue  of  stability  of 
predictive  control  designs  when  the  plant  model  is  uncertain,  such  as  the  work  of  Zafiriou 
(1990)  and  Genceli  and  Nikolaou  (1993).  There  is  an  interest  in  the  research  community  to 
revisit  the  predictive  control  design  techniques  with  the  intention  of  including  robustness 
features  that  guarantee  stability  or  adequate  performance  when  the  plant  model  is  uncertain. 
One  interesting  example  is  the  robust  quadratic  DMC  design  including  hard  constraints 
studied  by  Zafiriou  (1990).  This  work  uses  a contraction  mapping  first  proposed  by 
Economou  (1985)  to  derive  time-domain  conditions  for  robust  stability  with  respect  to 
uncertainty  in  the  impulse  response  coefficients  of  the  nominal  model,  but  this  approach 
involves  a very  large  numerical  computation  effort.  Genceli  and  Nikolaou  (1993)  propose 
an  analysis  and  synthesis  method  for  predictive  controllers  based  on  FIR  models,  including 
constraints  and  using  a linear  cost  functional.  These  authors  use  a parametric  model 
uncertainty  description  that  bounds  the  maximum  deviations  allowed  for  each  pulse- 
response  coefficient,  and  obtain  a sufficient  condition  for  robust  closed  loop  stability. 

The  robustness  of  predictive  controllers  designed  using  transfer  function 
representations  is  receiving  increasing  attention  in  the  literature.  Kouvaritakis  et  al.  (1992) 
propose  an  alternative  approach  to  GPC  that  employs  a precompensator  to  stabilize  the 
plant  before  the  predictive  design  is  carried  out.  The  (^-parameterization  procedure 
popularized  by  Youla  (Youla  et  al.  1976)  is  employed  in  order  to  construct  a final 
controller  that  is  robust  with  respect  to  unstructured  perturbations.  The  authors  state 
rigorous  necessary  and  sufficient  conditions  for  robust  stability;  however,  the  approach 
proposed  for  synthesizing  robust  controllers  is  an  approximate  albeit  practical  scheme.  The 
method  consists  of  using  polynomial  or  fixed-order  transfer  function  approximations  for 
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the  Youla  parameter,  and  least-squares  methods  to  identify  the  parameters  of  the  robust 
design.  Hrissagis  et  al.  (1996)  present  a direct  method  for  designing  a predictive  controller 
that  is  robust  with  respect  to  unstructured  perturbations.  In  this  approach,  an  appropriate 
Youla  parameter  is  explicitly  computed  by  solving  a model-matching  problem. 

This  chapter  presents  a method  for  designing  controllers  that  are  robust  with  respect  to 
real  parametric  uncertainty.  The  technique  is  based  on  the  results  of  Rantzer  and  Megretski 
(1994),  and  therefore  relies  on  using  fixed-order  approximations  for  the  Youla  parameter  in 
the  controller  parameterization.  The  nominal  predictive  controller  is  designed  using  well- 
established  methods  and  the  nominal  servo  performance  is  retained  by  the  robust  controller. 
Furthermore,  the  design  technique  can  be  easily  modified  to  incorporate  integral  action  in 
the  robust  controller,  allowing  for  the  rejection  of  asymptotically  constant  disturbances. 

The  chapter  is  organized  in  the  following  manner.  The  next  section  discusses  the 
design  of  a predictive  controller  of  the  GPC  type  for  a nominal  transfer  function  plant 
model.  The  third  section  details  the  parameterization  of  all  nominally  stabilizing  controllers 
through  the  use  of  the  Youla  parameter  Q.  The  robust  predictive  control  design  based  on 
the  ellipsoidal  uncertainty  description  is  discussed  in  the  fourth  section.  The  fifth  section 
contains  the  modifications  of  the  design  procedure  required  to  incorporate  integral  action 
into  the  robust  controller.  A further  modification  to  allow  the  inclusion  of  unstructured 
uncertainty  is  discussed  in  the  sixth  section.  An  example  is  given  in  the  seventh  section  to 
illustrate  the  proposed  design  method.  The  design  equations  used  to  construct  the  nominal 
predictive  controller  are  given  in  the  final  section. 

4.2  Nominal  Predictive  Control  Design 

Predictive  control  is  a control  methodology  that  is  well  documented  in  the  literature.  In 
particular,  a wealth  of  knowledge  is  available  to  resolve  crucial  design  issues  such  as 
nominal  closed-loop  stability,  and  parameter  tuning  (Lambert  1987;  Mohtadi  1987). 
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Predictive  controllers  are  usually  implemeted  by  executing  at  every  sampling  instant  an 
algorithm  that  solves  a quadratic  optimization  problem.  For  analysis  purposes,  it  is 
desirable  to  represent  the  algorithmic  controller  in  terms  of  transfer  functions,  allowing  the 
utilization  of  classical  z-domain  tools  for  analyzing  stability  and  performance.  This  section 
presents  a brief  review  of  the  analysis  technique  discussed  by  Crisalle  et  al.  (1989),  which 
casts  an  algorithmic  predictive  control  law  of  the  GPC  type  into  a form  involving  transfer 
function  operators.  The  resulting  nominal  controller  is  used  as  the  basis  for  the  design  of  a 
robust  controller. 

Consider  the  nominal  process  model 

A(z)  y(z)  = B(z)  w(z) 

where  y(z)  and  u(z ) are  the  process  output  and  input,  respectively,  and  A(z)  and  B(z)  are 
the  coprime  polynomials 


A(z)  = zn  + an-\Zn  *+•••  + 0o 

(4.2) 

B(z)  = bmz'n  + bm_\Zm  ’+...  + b0 

(4.3) 

of  order  n and  m,  respectively,  where  n>m.  Predictive  control  involves  the  selection  of 

future  control  moves  that  minimize  the  quadratic  cost  functional 

Ny  nu 

t)  = £[r(t  + i)  - y(t  + i|t)f  + X X[A«(t  + i)]2  (4.4) 

i=l  i=0 

where  {r(t+/)}  is  the  sequence  of  future  values  of  the  set  point,  {y(t+/lt)}  is  the  sequence 
of  predicted  future  values  of  the  output,  {Au(t+i)}  is  the  sequence  of  future  control 
increments,  A is  the  move-suppression  parameter  used  to  penalize  excessive  control 
energy,  and  parameters  Ny  and  Nu  are  the  prediction  and  control  horizons,  respectively. 
The  optimal  control  move,  the  w(t)  that  minimizes  the  functional  J(t)  for  the  prescribed  set 
point  sequence  { r(t) } , is  found  by  differentiating  (4.4)  with  respect  to  the  control  moves, 
equating  the  result  to  zero,  and  solving  for  «(t).  The  predictive  control  law  can  also  be  cast 
in  terms  of  transfer  function  operators  as  (Crisalle  et  al.  1989) 
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R(z)  / . tv  w x S(z) 
n u(z)  = T(z)r(z)  y(z ) 

z z 

(4.5) 

where 

R(z)  = zn  + rn-\Zn~X  +...  + r0 

(4.6) 

S(z)  = Snzn  + sn_i  z,!_1  + ...  + s0 

(4.7) 

x Nv  Nv- 1 

T (z)  = tN'fz  ' +tN y_jz  • +...  + qz 

(4.8) 

These  operators  satisfy 

ii 

o 

(4.9) 

and 

T(l)  = S(  1) 

(4.10) 

The  coefficients  of  the  moving-average  polynomial  S(z),  the  regressor  polynomial  R(z), 
and  the  set-point  advancement  polynomial  T(z)  are  functions  of  the  tuning  parameters  Ny, 
Nu,  and  A,  and  of  the  model  polynomials  A(z)  and  B(z).  Note  that  the  predictive  control 
law  (4.5)  includes  an  integrator  if  equation  (4.9)  is  satisfied.  A block-diagram 
representation  of  the  predictive  control  structure  is  shown  in  Figure  4.1a.  Specific 
equations  for  the  polynomials  (4.6)  - (4.8)  are  given  in  Section  4.8;  further  details  of  the 
derivation  can  be  found  in  Crisalle  et  al.  (1989).  A formulation  equivalent  to  (4.5)  is  also 
derived  in  McIntosh  et  al.  (1991). 

Note  that  the  transfer  functions  operating  on  u(z ) and  y(z)  in  the  nominal  predictive 
controller  (4.5)  are  biproper  and  of  order  n,  the  order  of  the  nominal  plant  model.  Note 
also  that  the  set-point  advancement  polynomial  T(z)  is  of  degree  equal  to  the  prediction 
horizon  Ny.  Since  Ny  > n is  a common  tuning  prescription  (Clarke  et  al.  1987),  the  order 
of  T(z)  may  exceed  the  order  of  R(z),  making  the  control  law  noncausal  with  respect  to  the 
set-point  signal.  This  noncausality  is  a natural  consequence  of  the  inclusion  of  future 
values  of  the  set  point  in  (4.4).  Figure  4.1a  shows  that  T(z)  acts  on  the  set  point  to 
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produce  the  intermediate  signal  v(z)=T(z)r(z),  which  has  the  simple  time-domain 
representation 

v(t)  = tNy  r(t  + Ny)  + tN^r(t  + Ny  - 1)  + ...  + hKt  + 1)  (4. 1 1) 

It  is  useful  to  remark  that  the  nominal  model  (4.1)  and  the  functional  (4.4)  are  simpler 
versions  of  more  elaborate  formulations  that  improve  the  design  performance  at  the  expense 
of  added  complexity.  Enhancements  of  the  predicitve  control  law  presented  above,  such  as 
the  inclusion  of  a lower  prediction  horizon  parameter  (Clarke  et  al.  1987),  the  addition  of  a 
weighted  end-point  term  in  (4.4)  to  guarantee  stability  for  arbitrary  parameter  choices 
(Kwon  and  Byun  1989;  Demircioglu  and  Clarke  1993),  and  the  use  of  a filtered  set  point, 
can  be  accommodated  within  the  framework  proposed  in  this  chapter  through  obvious 
modifications. 

Figure  4.1a  illustrates  the  closed  loop  system  established  when  the  nominal  predictive 
controller  (4.5)  is  connected  to  the  process  (4.1).  In  addition  to  the  set  point  signal  r(t),  the 
figure  also  shows  an  additive  output  disturbance  signal  d{ t).  Note  that  the  servo  dynamics 
of  the  closed  loop  are  fully  characterized  by  the  equations 

[ A(z)R(z)  + B(z)S(z)]  y(z ) = z"B(z)T(z)  r(z)  (4. 1 2) 

[A(z)R(z)  + B(z)S(z)]  u(z)  = z"  A(z)T(z)  r(z)  (4.13) 

Therefore,  the  stability  of  the  closed  loop  for  a given  nominal  predictive  controller  is 
contingent  on  the  location  of  the  roots  of  the  characteristic  polynomial  A(z)R(z)  + B(z)S(z). 
Furthermore,  due  to  the  presence  of  the  integral  action  (4.9)  in  the  controller  and  to  the  gain 
equality  (4.10),  the  closed  loop  dynamics  described  by  (4.12)  are  guaranteed  to  realize  zero 
offset  in  the  servo  response.  The  integrator  also  guarantees  perfect  steady-state  disturbance 
rejection  for  all  disturbance  signals  that  reach  a constant  steady-state.  These  desirable 
performance  characteristics  of  the  nominal  controller  will  be  preserved  in  the  robust 
predictive  controller  designed  in  the  following  sections. 
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Figure  4.1.  (a)  Structure  of  a nominal  predictive  controller,  (b)  Structure 
of  the  parameterized  predictive  controller  including  the  Youla  parameter  Q(z). 


4.3  Nominally  Stabilizing  Controller  Parameterization 

In  this  section,  all  nominally  stabilizing  controllers  are  parameterized  in  terms  of  the 
nominal  predictive  controller  (4.5)  and  a transfer  function  Q(z ) selected  in  the  spirit  of 
Wiener-Hopf  design  (Youla  et  al.  1976).  However,  a modification  in  the  parameterization 
is  introduced  to  achieve  two  important  design  requirements:  (i)  the  parameterized  controller 
must  preserve  the  servo  performance  and  the  steady-state  disturbance  rejection  properties  of 
the  nominal  controller,  and  (ii)  the  parameterized  controller  must  also  be  a predictive 
controller. 

Consider  a nominal  predictive  controller  (4.5)  that  stabilizes  the  closed  loop  system 
(4.12)  - (4.13).  Provided  the  closed  loop  is  stable,  the  nominal  closed-loop  characteristic 
polynomial 


A*(z)  = A(z)R(z)  + B(z)S(z) 


(4.14) 
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is  a Schur  polynomial  of  degree  2n.  As  a first  step  in  the  parameterization  of  all  nominally 
stabilizing  controllers,  a coprime  fractional  representation  of  the  nominal  plant  model  (4.1) 
is  constructed 


p0(t)  = B(z)=  Nizl 
y A(z)  M(z) 


(4.15) 


where  N(z)  and  M(z)  are  proper  and  stable  transfer  functions  that  satisfy  the  Diophantine 
equation 

A(z)X(z)+M(z)F(z)  = l (4.16) 


for  some  pair  of  stable  and  proper  transfer  functions  X(z)  and  F(z).  (Note  the  use  of 
italicized  capital  letters  for  transfer  functions,  while  polynomials  are  designated  with  plain 
capital  letters.)  A suitable  (M(z),  (V(z))  pair  is  readily  derived  from  the  nominal 
characteristic  polynomial  (4.14)  as  in  Hrisaggis  et  al.  (1996).  First,  the  closed-loop 
characteristic  polynomial  is  factored  into  the  form  A*(z)  = Ai(z)A2(z),  where  both  Ai(z) 
and  A2(z)  are  Schur  polynomials  of  degree  n.  If  A*(z)  contains  complex  poles  then  Ai(z) 
and  A2(z)  are  constructed  such  that  each  complex-conjugate  pair  is  contained  in  either  Ai(z) 
or  A2(z)  to  ensure  that  each  polynomial  factor  has  only  real  coefficients.  Both  sides  of 

(4.14)  are  then  divided  by  the  factored  characteristic  polynomial  to  obtain 

A(z)R(z)  | B(z)S(z)  (4.17) 

Aj(z)A2(z)  Aj(z)A2(z) 

Finally,  stable  and  proper  factorizations  satisfying  (4.16)  are  defined  as 


M(z)  := 


A(z) 
Aj(z)  ’ 


N(z):= 


B(z) 

Aj(z) 


(4.18) 


and 


nzy.=  RU) 


(4.19) 


A2(z)  A2(z) 

This  result  allows  the  nominal  predictive  control  law  (4.5)  to  be  written  in  the  equivalent 
form 


F(z)  u{z)  = Z(z)  r(z)  - X(z)  y(z) 


(4.20) 


where 
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(4.21) 


The  set  of  all  solutions  of  (4.16)  can  be  written  in  terms  of  the  transfer  functions 
(4.18)  - (4.19)  and  a stable,  proper  transfer-function  Q(z ) through  the  well-known  relations 
(Youla  etal.  1976) 


Therefore,  the  set  of  all  stabilizing  controllers  with  the  structure  (4.20)  is  parameterized  in 
the  form 


to  yield  the  control  structure  shown  in  Figure  4.1b.  Setting  Q(z)= 0 recovers  the  nominal 
predictive  controller  (4.20). 

In  contrast  to  the  standard  Youla  parameterization  approach,  the  transfer  function  X(z) 
+ M(z)Q{z)  appears  in  the  feedback  path  of  Figure  4.  lb,  instead  of  appearing  in  the  control 
block  immediately  preceding  the  plant.  This  approach,  adopted  from  Hrissagis  et  al. 
(1996),  in  conjunction  with  factorizations  (4.18)  and  (4.19)  that  make  use  of  the  nominal 
closed-loop  polynomial,  introduces  the  following  highly  desirable  properties  in  the 
parameterized  input-output  maps. 

Property  4.1.  The  nominal  control  loop  of  Figure  4.1a  and  the  parameterized 
control  loop  of  Figure  4.1b  have  identical  servo  transfer  functions  y(z)/r(z)  and  u(z)/r(z). 

Proof.  The  proposition  is  proved  by  carrying  out  block-diagram  algebra  on  each 
figure  to  derive  in  both  cases  the  servo  transfer  functions  y(z)/r(z)  and  u(z)/rfz).  ♦ 

Property  4.2.  Given  that  the  nominal  controller  (4.5)  is  a predictive  controller,  then 
the  parameterized  controller  (4.24)  is  also  a predictive  controller. 


X\z)  = X(z)  + M(z)Q(z) 


(4.22) 


nz)=Y(z)-N(zMz) 


(4.23) 


[Y(z)  - N(z)Q(z)]  u(z)  = Z(z)  r(z)  - [X(z)  + M(z)Q(z)  ] y(z)  (4.24) 


Proof.  If  (4.5)  is  a predictive  controller,  then  by  definition  it  yields  a control  sequence 
{ w(t) } that  minimizes  the  predictive  performance  index  (4.4)  for  any  prescribed  set-point 
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trajectory  { r(t) } . For  the  given  set  point  trajectory,  it  follows  from  Property  4.1  that  the 
parameterized  controller  (4.24)  will  also  produce  the  same  control  sequence  due  to  the 
equality  of  the  servo  transfer  function  u(z)/r(z).  It  follows  that  the  parameterized  controller 
is  also  a predictive  controller  because  it  yields  a control  sequence  that  minimizes  (4.4).  ♦ 

Since  any  allowable  parameter  Q(z ) yields  the  same  servo  transfer  functions  y(z)/r(z) 
and  u(z)/r(z),  the  parameterized  controller  has  the  intrinsic  capability  of  preserving  the 
nominal  servo  performance.  Also  note  that  although  the  terms  containing  Q(z)  effectively 
cancel  out  in  the  servo  transfer  functions,  the  transfer  function  e(z)/v\(z)  = M(z)[Y(z)  - 
N(z)Q(z)]  in  Figure  4.1b  is  affine  in  Q(z),  as  in  the  standard  Youla  parameterization 
method. 


4.4  Robust  Predictive  Control  Design 


In  order  to  incorporate  uncertainty  into  the  plant  description,  the  nominal  plant  transfer 
function  shown  in  Figure  4.1b  is  now  represented  as  the  uncertain  transfer  function 
B(z)  + AB(z)  _ bmz'n  + ...  + b0 +bbmzm  + ...  + bbo 


P{Z ) = 


(4.25) 


A(z)  + AA(z)  zn  + an-\Zn  +...  + #o  +ban_\Zn  +...  + 5ao 
where  A(z)  and  B(z)  are  given  in  (4.2)  and  (4.3)  and 

AB(z)  = bbmzm  + ...  + bbQ  (4.26) 

AA(z)  = 8an_1zn_1+...  + 5fl0  (4.27) 


The  values  of  the  coefficients  of  AA(z)  and  AB(z)  are  not  known  explicitly,  however,  the 
perturbation  vector 

8p  = [5a„_1  ...  daQbbm  ...  8/?0]T  e cr'!+"!+1  (4.28) 


composed  of  the  coefficients  of  AA(z)  and  AB(z),  is  constrained  to  lie  in  an  ellipsoid  Up  in 
the  parameter  space.  This  ellipsoid  is  defined  by  a positive  definite,  symmetric  matrix  Qp 
such  that 
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Ep  = j8p  e 9T+m+1|  5pTQ”18p  < l} 


(4.29) 


The  objective  of  the  robust  control  design  is  to  select  a suitable  Youla  parameter  Q(z)  such 
that  the  resulting  controller  stabilizes  the  family  of  plants  described  by  Ep. 

The  work  of  Rantzer  and  Megretski  (1994)  considers  robust  controller  design  for  any 
system  where  the  uncertainties  can  be  extracted  from  the  closed  loop  as  shown  in  Figure 
4.2,  where  it  is  noted  that  the  signal  vv  is  scalar.  This  formulation  immediately  restricts  this 
technique  to  rank  one  M-A  problems.  The  affine  uncertainty  description  in  equation  (4.25) 
results  in  a rank  one  M-A  problem,  so  this  formulation  can  be  accommodated  by  the 
technique  proposed  by  Rantzer  and  Megretski  (1994).  It  is  straightforward  to  transform 
the  system  shown  in  Figure  4.1b,  now  with  an  uncertain  plant,  into  the  form  of  Figure  4.2. 
This  is  done  by  first  extracting  the  uncertainty  vector  8p  out  of  the  plant  block,  and  then 
constructing  the  appropriate  closed  loop  transfer  function  Tzw  from  the  uncertainty  block 
output  w to  the  uncertainty  block  input  z.  The  first  step  in  this  process  is  achieved  by  first 
rewriting  the  plant  block  as  shown  in  Figure  4.3. 


w 


Figure  4.2.  Feedback  structure  from 
Rantzer  and  Megretski. 


w 


5pT 


Figure  4.3.  Augmented  plant  with 
uncertainty  vector  8p. 


For  the  model  given  in  equation  (4.25),  Figure  4.4  shows  explicitly  the  internal  structure  of 
the  augmented  plant  given  in  Figure  4.3.  The  closed  loop  shown  in  Figure  4.3  is  given  by 


z 

'Gil  U) 

Gi2(z)" 

w 

_y_ 

_G2i(z) 

G22(z)_ 

u 

(4.30) 
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From  Figure  4.4,  it  is  straightforward  to  show  that  both  Gn(z)  and  G\2(z ) are  transfer 
function  vectors  with  n+m+1  rows  while  G2i(z)  and  G22U)  are  scalar  transfer  functions. 
The  exact  form  for  these  transfer  functions  is  given  below. 


Gn(z) 


A (z) 
~1 

A(z) 

0 


0 


Gji(z)  = -1 


G\2(z) 


z/,~1B(z) 

A2(z) 

B(z) 
A2(z) 
zm 
A (z) 


G22(z) 


1 

_ ~A(z) 

B(z) 

A(z) 


(4.31a) 


(4.31b) 


Figure  4.4.  Augmented  plant  in  detail. 

When  the  plant  block  in  Figure  4.1b  is  replaced  by  that  shown  in  Figure  4.3,  the  resulting 
closed  loop  is  transformed  to  the  form  shown  in  Figure  4.5.  The  overall  closed  loop  from 
w to  z becomes  Tzw,  where  the  transfer  function  vector  Tzw  is  given  by 
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TIW  = 7j  + T2Q  = {G„-G,2MXG21}  + {-G,2M2G2,}e  (4.32) 

The  expressions  in  (4.32)  for  the  transfer  function  vectors  T\(z ) and  72(z)  simplify  to 


z"_1R(z)  ' 

A1(z)A2(z) 

A?(z) 

R (z) 

B (z) 

A,(z)A2(z) 

and  T2(z)  = 

A?(z) 

zmS(z) 

zmA(z) 

Aj(z)A2(z) 

A\  (z) 

s'(z) 

A (z) 

Ai(z)A2(z)_ 

A f(z)  _ 

Note  that  all  elements  of  both  7i(z)  and  T2(z)  are  stable,  strictly  proper  transfer  functions, 
and  that  both  Ti(z)  and  Tjiz)  have  n+m+1  rows. 


Figure  4.5.  RST  configuration  with  the  augmented  plant. 


w 


z 


Figure  4.6.  Closed  loop  from  w to  z. 


65 


In  the  controller  design  method  of  Rantzer  and  Megretski  (1994),  the  approach  is  to 
express  the  Youla  parameter  Q(z)  as  a ratio  of  stable  proper  transfer  functions  (3(z)  and  a(z) 
in  the  form 

Q{z)  = P(z)/a(z)  (4.34) 

The  main  result  of  Rantzer  and  Megretski  (1994)  is  that  the  closed  loop  in  Figure  4.6  is 
stable  for  perturbations  5p  satisfying  I8pl  < 1 if  a(z)  and  (3(z)  satisfy 

|Re{7j (ej0} )a(eJ(0 ) + T2 (e;co )P(e7“ )}|^  < Re{a(<?yco)}  (4.35) 

for  all  to  e [0, 7i] . The  norm  used  in  measuring  the  size  of  5p  is  any  appropriate  vector 
norm,  and  |»|^  is  the  dual  norm  to  this  norm  (Luenberger,  1969). 

Since  the  nominal  servo  performance  of  the  predictive  controller  is  preserved  for  any 
choice  of  Q(z),  a possible  method  for  choosing  this  parameter  would  be  to  maximize  the 
stability  robustness  of  the  system  with  respect  to  the  size  of  the  uncertainty  vector  8p. 
Note  that  if 

|Re{r1(ey®)a(^“)  + r2(ey“)P(^C0)}|'/  < Reject®)}  for  |5p|<l  (4.36) 

where  y > 1 , then 

|Re{r1(^“)a(^'w)  + r2(e7'C0)p(^“)}|J  < Re{a(e7'“)}  for  |5p|  < y (4.37) 


Therefore,  if  a (z)  and  |3(z)  are  selected  to  solve  the  optimization 

Re{7](e7CO)a(ey“)  + r2(e7'(D)p(e7'(0)}|< 


min  max  <( 

cc.peRH^,  ooe[0,7t] 


Re 


{a(<?;co)} 


(4.38) 


subject  to  the  constraint 

Reja(<?7C0)j  >0  Vco  e [0, 7t] 


(4.39) 


then  the  resulting  closed  loop  will  be  robustly  stable  for  the  largest  possible  value  of  I8pl. 
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Since  a(z)  and  (3(z)  can  be  of  any  order  (such  that  the  resulting  Q(z)  is  proper),  this  is 

an  infinite  dimensional  optimization.  This  optimization  is  converted  to  a finite  dimensional 

problem  by  restricting  a(z)  and  (3(z)  to  be  of  the  form 

N 2N+1 

oc(x,z)  = 1 + %xkak(z)  and  (3(x,z)  = X'TtMz)  (4.40) 

*=i  k= N+l 


where  x = [^  •••  jc2n+i]T  e ^2N+I- 

Quiz) 


Then  Q(z)  is  a rational  transfer  function  of  order  N 
2N+1 

:=  (4.41) 

i + 

k= i 


The  functions  CLk(z)  and  Pj.(z)  should  be  chosen  such  that  the  limiting  (?N(z) 

lim  0N(z)  = lim  {(3(x,z)/a(x,z)} 

N— »<x>  N— 


(4.42) 


can  approximate  any  stable  transfer  function  Q(z ) arbitrarily  well.  Here,  the  basis  functions 
a k(z)  and  P^(z)  are  chosen  as 


ak(z)  = P*(z) 


( -az  + 1 ") 


k 


(4.43) 


V z-a  y 

where  a is  taken  to  be  a real  scalar  that  satisfies  Id  < 1.  This  choice  constrains  the  basis 
functions  to  be  all-pass  functions.  An  interesting  choice  for  a is  simply  a = 0,  reducing  the 
basis  functions  to  the  form 


a kiz)  = Mz)  = z 


(4.44) 


Obviously,  a ratio  of  polynomials  in  z_1  can  recover  any  rational  transfer  function  as  N 
grows  without  bound. 

Choosing  basis  functions  of  the  form  (4.43),  the  functions  a(x,z)  and  (3(x,z)  become 


N 


a(x,z)  = 1 + Yjxk 

k= 1 


(~az  + 1 


\k 


\ z-a  y 


2N+1 

and  P(x,z)  = X xk 


( -az  + 1 


lt-N+1 


(4.45) 


*=n+i  v z-a  y 

Then,  the  closed  loop  in  Figure  4.6  is  stable  for  all  allowable  perturbations  8p  (that  is,  all 
8p  satisfying  I8pl  < 1)  if  a(x,z)  and  P(x,z)  satisfy 
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Re{r1(^“)a(x,^w)+r2(^“)p(x,eyw)}|'/  < Re{a(x,eya))}  (4.46) 

The  constraint  in  equation  (4.29)  is  SpTQ“’8p  < 1.  A new  vector  8 can  be  defined 

5 :=  Qp1/2Sp  (4.47) 

i j j 

where  the  symmetric  positive  definite  matrix  Qp  satisfies 

Q,  = Qp2Qp'2  (4.48) 

Then  the  vector  8 satisfies  ||8||.,  < 1 for  all  allowable  8p.  The  block  diagram  shown  in 
Figure  4.3  can  be  adjusted  to  use  the  new  vector  8 as  the  perturbation  vector,  but  a factor 
Qj/2  must  be  included  in  the  plant  matrix  G(z),  as  shown  in  Figure  4.7.  The  new  vector  z 
is  related  to  z simply  by  z = Qp  z.  Therefore,  the  transfer  vectors  G\\  and  G \2  are 
changed  to  Qj,/2Gn  and  Qj,/2Gi2,  respectively.  Plugging  these  new  expressions  into 
equation  (4.31)  yields 

Tiw  = Q 'i\w  = 27j  - Q^'27-2e  (4.49) 

Note  that  all  the  elements  of  the  matrix  Q^,/2  are  real,  and  that  the  dual  norm  to  the 
Euclidean  2-norm  is  the  Euclidean  2-norm  itself.  Therefore,  equation  (4.35)  becomes 
|q|/2  Re|r1(e7“)a(x,e7W)-i-  r2(e-/<0)|3(x,e7“)|||2  < Re{a(x,r*“)}  Va)e[0,7t]  (4.50) 


Figure  4.7.  Augmented  plant  with  new  uncertainty  vector  8. 


Alternatively,  the  norm  used  to  measure  8p  can  be  defined  as 


|6p|  :=  q:i,2sp 


and  then  the  dual  norm  in  equation  (4.35)  is 


Ixl*  :=' 


q!/2x 


(4.51) 


2 


(4.52) 
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which  yields  exactly  the  result  shown  in  equation  (4.50). 

Equation  (4.50)  cannot  be  satisfied  if  Reja(x,eyt0)j  < 0 for  any  frequency. 
Therefore,  provided  that 

min  Re{a(x,e7a))}  > 0 (4.53) 

coe[0,7t]  *•  J 


then  the  closed  loop  is  robustly  stable  (for  a specific  value  of  x)  if 


max 

(flg[0,7t] 


Qj,/2  Re{7i(e^)a(x,^'®)  + T2(ej(°Mx,ej<0)} 


Re|a(x,e7<0)| 


< 1 


(4.54) 


As  outlined  above,  then  the  optimally  robust  controller  (actually  the  optimal  vector  x)  can 
be  found  by  solving  the  optimization  program 


xopt  = arg  min  max 
coe[0,7t] 


q;/2  r e{r,(Oa(xyw)  + r2(^)P(x,^)}|| 


Re 


{a(x,e7C0)} 


(4.55) 


subject  to 


min  Re{a(x,e-yC0)}  > 0 
CO€[0,7t]  ^ J 


(4.56) 


for  larger  and  larger  values  of  N.  In  practice,  a value  of  N will  be  found  such  that  the  value 

Qp2  Re{r1(e;“)a(x,e;“)  + r2(e-/0))P(x,e-/0))}| 


max 

(O6[0,7t] 


Re 


joc(x,e;to)j 


(4.57) 


does  not  decrease  appreciably  as  N is  increased  further.  The  optimal  vector  x°Pl  for  this 
value  of  N is  considered  to  describe  the  optimal  controller.  Alternatively,  the  value  of  N 
can  be  chosen  to  enforce  a maximum  degree  condition  on  the  resulting  controller.  The  final 
form  of  the  robust  predictive  controller  is 


CN(z) 


X(z)  + M(z)<2n(z) 

Y(z)-  N(z)Qti(z) 


(4.58) 
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which  will  be  biproper  and  of  degree  2n+N  if  there  are  no  common  factors  that  cancel  in 
forming  Cn(z).  Thus,  a specific  value  for  N can  be  chosen,  and  the  "optimal"  QN(z) 
found  for  that  value  of  N. 

The  use  of  basis  functions  of  the  form  (4.45)  results  in  the  following  desirable 
properties  of  the  optimization  (4.55)  - (4.56)  as  outlined  below. 

Lemma  4.1.  If  the  basis  functions  a(x,z)  and  (3(x,z)  are  given  by  equation  (4.45), 
then  the  objective  function 


£(x)  :=  max 

coe[0,7t] 


Qp2  Re{7j  (ej(0  )a(x,  ej®  ) + T2  (ej(°  )P(x,  ej(0 )} 
Reja(x,e7C0)j 


(4.59) 


is  a quasi-convex  function  of  x,  and  the  constraint  function 


0(x)  :=  min  Re 

a)e[0,7t] 


{a(x,<?7“)} 


(4.60) 


is  a concave  function  of  x. 

Proof.  The  proof  is  given  in  the  Appendix. 


♦ 


Any  general  optimization  routine  for  quasi-convex  functions  requires  gradients  or 
subgradients  of  both  the  objective  and  constraint  functions.  It  is  noted  that  for  any  given 
frequency  it  is  possible  to  write 


Qj/2R  e{  r,  (e^  )cx(x,  ) + r2  (^^“  )P(  x,  )}||2  ||q^/2{A(o))x  + b(co)}| 


Re 


{a(xy“)} 


1 + c (oo)x 


(4.61) 


where  the  matrix  A(co)  and  vector  b(co)  depend  on  7](e7CO),  T2(ejUi),  and  the  pole  of  basis 
function,  and  vector  c(co)  depends  on  the  pole  of  basis  function.  If  the  value  to  =to*  that 
satisfies 


max 

coe[0,7i] 


Q^/2Re{7](e7“)a(x,e7“)  + r2(^“)P(x,e7(0)}[ 

Rejalx,^®)! 


Qp2{A(0)*)x  + b((0*)}||2 
l + cT(co*)x 


(4.62) 


is  found,  then  an  appropriate  subgradient  (x)for  the  objective  function  (4.59)  is 
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d [ Qp2{A(co*)x  + b(a>*)} 
dx  1 + cT(co*)x 


*(x) 


(Boyd  and  Barratt,  1991).  Similarly,  if  the  value  to  =£0*  that  satisfies 

min  Reja(x,e;w))  = 1 + cT(co*)x 

roe[0,7t]  *•  J 


(4.63) 


(4.64) 


is  found,  then  an  appropriate  subgradient  ^(x)  for  the  constraint  function  (4.60)  is 

*00  = ~(l  + cT(0)*)x)  (4.65) 

To  compute  the  subgradients  for  the  constraint  and  objective  functions,  the  extreme 
frequencies  that  satisfy  equations  (4.62)  and  (4.64)  must  be  found.  The  details  of  this 
process  are  discussed  in  the  following  subsection. 


4.4.1  Constraint  Testing  and  Objective  Function  Value  Computation 

The  optimization  (4.55)-(4.56)  can  be  performed  using  standard  methods  for 
constrained  quasi-convex  optimization,  such  as  the  ellipsoidal  method  outlined  in  Boyd  and 
Barratt  (1991).  This  method  requires  the  value  of  the  objective  function  and  its  subgradient 
if  the  current  point  is  feasible  (satisfies  the  constraint  (4.56)),  and  the  value  of  the 
constraint  and  its  subgradient  if  the  constraint  is  violated.  As  commented  above, 
computation  of  either  subgradient  requires  finding  the  frequency  that  satisfies  (4.62)  or 
(4.64).  The  case  of  the  objective  function  is  considered  first,  since  the  methods  used  are 
very  similar  to  the  methods  introduced  in  Chapter  3. 

If  the  constraint  (4.56)  is  satisfied,  then  the  function 

|Qj/2Re{r1(c^)a(x,c^)  + r2(c^)P(x,eJto)}||2 
UX>  :=  Re{a(x.O}  <4  66> 

is  non-negative  for  all  (0  e [0,7t].  Therefore,  it  is  possible  to  find  a stable  rational  transfer 
function  whose  magnitude  matches  the  value  of  ^(x)  for  all  co  e [0, 7c] , using  the  methods 
of  Chapter  3.  Then  the  objective  function  value  is  found  as  the  Hoo  norm  of  this  transfer 
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function,  and  it  is  also  possible  to  compute  the  maximizing  frequency  so  that  the 
subgradient  can  be  computed.  There  are  several  additional  points  that  must  be  addressed 
before  the  method  of  Chapter  3 is  used  to  find  the  transfer  function  fit  to  (4.66).  The  first 
is  that  only  the  real  part  of  T[(e-/“)a(x,e-/a))  + r2(e-/“)P(x,e-/0))  is  retained.  Since  the 
quantity  7’1(e'/CO)oc(x,e-/a))  + ^(e-^Plx,^03)  is  a transfer  function  vector,  this  is 
equivalent  to  computing  the  real  part  of  each  element.  Given  any  general  polynomial  in  z 

f(z)  = fnpz"P+-  + fo+-  + f-nnz'""  (4.67) 

then  the  real  part  of  this  polynomial,  evaluated  along  the  unit  circle,  is 

Rejfte^)!  = fn  cos(npco)  + ...+  f0  + .,.  + f_„n  cos(nnco)  (4.68) 


and  if  it  assumed  that  nn>np,  then 

Re  jf(e-/0))J  = fo  + (f_,  +f1)cos(o))  + ...  + (fnp  + f_np)cos(npco)- 

f(_n  l)Cos((np  + l)co)  + ...  + f_nn  cos(nnco) 


(4.69) 


This  is  equivalent  to  the  frequency  response  magnitude  of  the  polynomial 
fr(z):=-f_nnznn+...  + -f(  _ i)Z  p +-(f_„  +fn  )z  p +...  + -(f_i +fi)z  + f0 + 


-(f_l+fl)z  1+...  + |(f_n  +fn  )z  np+|f(-n  -\)Z  "P  1+...  + ^-f_nnZ  "n 


(4.70) 


where  it  is  noted  that  the  frequency  response  of  (4.70)  is  purely  real  by  definition. 
Therefore,  it  is  possible  to  compute  the  coefficients  of  a polynomial  whose  frequency 
response  magnitude  is  equivalent  to  the  real  part  of  any  individual  element  of  the  vector 
T\(ej(i>)a(x,e-i(i))  + 7’2(e-/W)P(x,e-/0)),  and  all  of  these  elements  can  be  stacked  to  forma 
matrix  of  coefficients.  The  effect  of  the  matrix  multiplication  by  Q^/2  is  accounted  for 
simply  by  multiplying  the  matrix  of  coefficients  on  the  left  by  Q|/2 . Each  polynomial  is 
then  multiplied  with  itself  to  achieve  the  effect  of  the  squaring  involved  in  the  two-norm,  all 
of  the  coefficients  corresponding  to  like  powers  are  summed,  and  finally  a spectral  factor  is 
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found  for  the  numerator  to  account  for  taking  the  square  root.  Before  this  is  done,  a 
common  denominator  for  the  two-norm  part  and  Re|a(x,e-/<0)|  is  found  and  extracted.  A 
spectral  factor  for  Reja(x,c-/C0)j  is  also  found,  but  it  is  multiplied  by  itself  since  there  is 
no  need  to  take  the  square  root  of  Re|a(x,cj,0))| . 

The  constraint  (4.56)  can  be  tested  using  the  Kalman- Yakubovich-Popov  lemma.  The 
form  of  this  lemma  given  in  Rantzer  (1993)  is  reproduced  here: 

Lemma  4.2.(Kalman-Yakubovich-Popov)  Given  G(s)  = C(sI-A)_1B  + D 
with  A stable,  the  following  statements  are  equivalent: 

(i)  G(s ) + G*(s ) < 0 forallRes>0 

(ii)  There  exists  a positive  definite,  symmetric  matrix  P such  that 


p 

o' 

A 

B 

At 

CT1 

p 

O' 

0 

I 

C 

D 

+ 

bt 

dtJ 

0 

i 

< 0 


(4.71) 


For  a scalar  G(s),  condition  (i)  becomes  Re(G(s))  < 0.  For  the  discrete  time  case,  the 
transfer  function  cx(x,z)  is  transformed  into  continuous  time  by  means  of  a bilinear 
transformation.  Define 

1 + 5^ 


a(x,z)  = a x, 


1 -s) 


:=  a'(x,s) 


(4.72) 


and  note  that 


a 


i-rn 

X, 


V 


1 + Z 


-1 


= a(x,z) 


If  Re{a'(x,y(o)}  > 0 for  to  > 0 implies  Re 


a 


x, 


l-e’*' 


1 + e 


-ye 


since 


7 


(4.73) 

> 0 for  0e[O,7t],  then. 


Re 


[ ( i-,-y0YI 


a 


x, 


1 + e 


-ye 


= Re 


M^e)} 


(4.74) 


Lemma  4.2  can  be  used  to  test  Re{a'(x,yco)}  > 0 instead  of  Re  ja(x,eyC0  j| . Now 


\-e  1 -h 2 / sin 0 — 1 / sin0 

= J 


1 + e 


-ye 


2 + 2cos0 


1 + COS0 


(4.75) 
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so  if  co:= 


( sin0 


, then  as  0 varies  from  0 to  n,  to  varies  from  0 to  Therefore,  with 


Vl  + cos0 

this  definition  of  co,  Reja'(xjco)}  > 0 for  to  > 0 implies  Re^,  va  ^ 

1 + e ^ j 

0 e [0,Tt].  The  transfer  function  a'(x,s ) is  a scalar,  but  we  wish  to  test  Re(cc'(x,s))  > 0 so 


a 


x, 


> > 0 for 


we  take  a'Cx,^  = -G(s).  A state  space  description  of  a'(x,s)  can  be  found  easily  in 
controller  canonical  form,  and  the  state  matrix  A will  be  stable  by  definition  of  the  basis 
functions  chosen.  Given  the  state  space  description  of  a\x,s ) as  [A,B,C,D], 
Re(  a'(x,s))  > 0 for  all  Re  s > 0 becomes  equivalent  to  the  existence  of  a positive  definite, 
symmetric  matrix  P such  that 


p 

o' 

A 

B 

AT 

-CT1 

p 

o' 

0 

I 

-C 

-D 

+ 

bt 

-DtJ 

0 

i 

< 0 


(4.76) 


This  condition  is  a linear  matrix  inequality  (LMI)  feasibility  problem  in  the  matrix  variable 
P,  and  (4.76)  can  be  checked  using  the  LMI  toolbox  in  Matlab.  This  formulation  of  the 
constraint  is  very  useful  in  the  continuous  time  version  of  this  controller  design  method. 


4.5  Robust  Control  Design  With  Steady  State  Disturbance  Rejection 

Offset-free  regulation  in  the  presence  of  asymptotically  constant  disturbances  will  not 
be  attained  in  general  from  the  robust  controllers  designed  using  the  methods  of  the 
previous  section.  In  this  section,  a simple  modification  derived  in  Hrissagis  et  al.  (1996) 
is  discussed  that  ensures  that  the  robust  controller  obtained  provides  integral  action.  The 
nominal  predictive  controller  (4.5)  leads  to  the  nominal  regulation  transfer  function 

y(z)  = A(z)R(z) 

d(z ) A(z)R(z)  + B(z)S(z)  1 ' 

from  which  it  follows  that  limy(f)  = 0 for  step  disturbances  as  well  as  for  other 

/— >00 

disturbances  with  a constant  steady  state  because  R(  1)  = 0.  On  the  other  hand,  the  nominal 
regulation  transfer  function  for  the  parameterized  robust  controller  (4.24)  is 

y(z)  _ A(z)R(z)  A(z)B(z) 


d(z)  A(z)R(z)  + B(z)S(z)  A?(z) 


■Q(z) 


(4.78) 
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Because  the  synthesis  procedures  described  in  the  previous  section  do  not  necessarily  yield 
a Youla  parameter  satisfying  Q{\)  = 0,  the  robust  predictive  controller  may  display 
unacceptable  nominal  regulation  performance  at  the  steady  state,  unless  the  nominal  plant 
has  an  integrator  (that  is,  A(1)=0)  or  is  a self-regulating  process  (B(l)  = 0). 

Clearly,  the  robust  predictive  controller  will  attain  perfect  steady  state  disturbance 
rejection  for  all  the  plants  belonging  to  the  uncertain  family  (4.25)  only  if  the  Youla 
parameter  has  a zero  gain,  that  is,  2(1)  = 0.  This  gain  constraint  can  be  introduced  in  the 
robust  predictive  controller  design  through  a simple  modification  of  the  factorizations 
(4.18)  - (4.19).  First,  the  integrator  is  extracted  from  the  nominal  predictive  controller  by 

writing  R(z)  = (z-l)R'(z),  and  then  (4.17)  is  rewritten  in  the  form 

A(z)(z-l)R'(z)  , B(z)S(z) 


Aj(z)A2(z)  ' Aj(z)A2(z)  1 
Introducing  the  modified  coprime  factorization 

M(2):=Ozi)A«  ^BU) 


and 


X(z):= 


zAj(z) 

S(z) 


Y(z):= 


A(z) 
zR'(z) 


(4.79) 


(4.80) 


(4.81) 


A2(z)  A 2(z) 

A A A A 

leads  to  operators  that  satisfy  the  Diophantine  equation  N(z)X(z)+M(z)Y(z)= 1.  Note  that 

this  modified  factorization  is  equivalent  to  augmenting  the  nominal  plant  with  an  integrator. 

z 

Thus  the  modified  plant  output  is  y := y so  that  (4.30)  becomes 


z — 1 

Gii(z)  G\ 2(z) 


w 


y\  [g2\(z)  G22(z)]lu_ 

(4.82) 

where 

G2\(.z)  = — ^—G21(z) 
z — 1 

(4.83) 

G22(z)  = G22{z) 

z — 1 

(4.84) 

These  modified  transfer  functions  are  then  substituted  into  (4.32)  resulting  in 
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71  = G|  | — Gl2MXG2l  = G]]—G] 


T2  — — G^Af^G^l  — — G\ 


1 -'-712 


(iziW_s ' 

z J Vz-17 
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V Z 


z-\) 


G2i  - 


G21  = 71  (4.85) 

fz-lV 


\ Z ) 


(4.86) 


It  is  now  possible  to  solve  for  a parameter  £>N(z)  using  the  design  procedure  detailed  in  the 
previous  section.  After  a solution  is  found,  the  Youla  parameter  £?n(z)  used  in  the 
parameterized  predictive  control  structure  of  Figure  4.  lb  is  constructed  by  re-associating 
the  augmented-plant  integrator  with  the  controller  to  obtain 

(4.87) 


2n,i(z)  = ^^<2n(z) 


The  final  robust  predictive  controller  design  for  this  case  is  obtained  by  substituting  the 
Youla  parameter  (4.87)  and  the  factorizations  (4.18)  - (4.19)  into  the  structure  (4.24).  The 
resulting  controller  includes  an  integrator  since  (4.87)  satisfies  the  zero-gain  condition 
0n(1)  = O. 

4.6  Inclusion  of  Unstructured  Uncertainty 


High-frequency  dynamics  that  are  often  neglected  in  system  modeling  can  be 
represented  in  most  cases  as  unstructured  uncertainty.  In  this  section,  unstructured 
uncertainty  is  incorporated  into  the  uncertainty  description  and  a modified  algorithm  for 
controller  design  is  presented. 

The  structural  constraint  on  the  uncertainty  description  in  Figure  4.2  is  that  the  signal 
w must  be  a scalar.  Therefore  it  is  possible  to  include  an  element  A(z)  along  with  the 
vector  8p  in  the  uncertainty  block,  provided  w remains  a scalar.  Figure  4.8  shows  one 
possible  combination.  This  structure  is  similar  to  the  additive  uncertainty  model  that 
appears  in  the  literature  in  that  if  all  the  elements  of  8p  are  zero,  then  the  plant  block  is 
equivalent  to  the  transfer  function 


P(z)  = 


B(z) 

A(z) 


+ A(z) 


(4.88) 
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or 


P(z)  = P°(z ) + A(z) 


(4.89) 


Figure  4.8  Plant  block  with  unstructured  element  A(z)  included. 


This  is  the  standard  additive  uncertainty  model  provided  ||A(z)||oo  < 1.  Obviously,  an 
appropriate  stable,  proper  weighting  function  can  be  extracted  from  any  stable  A(z)  as 
A(z)  = W(z)A(z)  so  that  ||A(z)||oo  < 1. 

When  the  elements  of  8p  are  non-zero,  then  the  general  transfer  function  from  u to  y is 


P(z) 


B(z)  + AB(z)  + A (z) 

A(z)  + AA(z)  A(z)  + AA(z) 


(4.90) 


and  a suitable  weighting  function  can  be  extracted  such  that 


P(z) 


B(z)  + AB(z)  + 
A (z)  + AA(z) 


A(z) 

A(z)  + AA(z) 


W(z)A(z) 


(4.91) 


which  satisfies  ||A(z)|  <1. 

The  inclusion  of  the  unstructured  element  A(z)  changes  the  stability  condition  (4.35)  to 
the  following  form  (Rantzer  and  Megretski,  1994) 
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Qp2  Re{7i5(e-/(D)a(x,e-/W)  + 7^  (^■/C°)P(x,e-/C0)|j^ 


+ 


r1A(e;l0)a(x,^l0)  + r2A(e;'0))(3(x,^0))|  < RejaCx,^10)}  (4.92) 


£ ? 

for  to  g [0, 7t] . The  transfer  function  vectors  T°(z)  and  T2 (z)  are  exactly  the  Ti(z)  and 
Ti(z)  given  in  (4.33),  while  7jA(z)  and  T2A(z)  are  simply 

7|a(z)  = -W(z)M(z)X(z)  and  T2A(z)  = -W(z)M2(z)  (4.93) 


4.7  Examples 


4.7.1  Example  1 


Consider  the  stable  nominal  plant  model 

p(z)  = ~— — (4.94) 

Z - .4 

a m 

with  the  corresponding  nominal  plant  parameter  vector  p =[-.4  5]  . The  matrix  Qp 
defining  the  plant  parameter  uncertainty  ellipsoid  is 


-0.09 

1.0 


(4.95) 


The  allowable  uncertainty  region  described  by  this  ellipse  is  shown  in  Figure  4.9.  It  is 
noted  that  the  range  a^  g (-1,1)  defines  all  stable  plant  models. 

A nominal  predictive  controller  is  designed  for  the  system  using  the  tuning  parameters 
Ny  = 1,  Nu  = 1,  and  X = 0,  resulting  in  the  following  controller  and  prefilter  polynomials 

R(z)  = z-1  (4.96) 

S(z)  = 0.28z-0.08  (4.97) 


T(z)  = 0.20z 


(4.98) 
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Figure  4.9  Allowable  uncertainty  region  for  the  example. 


It  is  readily  verified  that  these  polynomials  satisfy  (4. 10)  - (4. 1 1).  The  nominal  closed 
loop  polynomial  is 

A *(z)  = z2  (4.99) 

and  the  polynomials  Aj(z)  and  A2(z)  are  taken  as  Ai(z)  = z and  A2(z)  = z.  For  this 
example,  the  pole  of  the  basis  function  is  taken  to  be  a = 0.  It  is  straightforward  to  show 
that  the  perturbation 

8p  = [0.5  0.25]t  (4.100) 

yields  5p  8p  = 0 .9978  and  is  therefore  allowable.  The  plant  model  resulting  from  this 
perturbation  is 

5 25 

P(z)=—~  (4.101) 

z +.1 

The  closed  loop  polynomial  resulting  from  this  plant  model  and  the  nominal  controller  is 

G(z)  = z2  + 0.57z  - 0.52 


(4.102) 
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which  has  a root  outside  the  unit  circle,  showing  that  the  nominal  controller  is  not  robustly 
stable. 

Robust  controllers  with  and  without  integral  action  are  designed  for  the  choice  N = 1. 
The  optimal  vector  x for  the  robust  control  design  without  integral  action  is 

xopt  = [-0.0008  -0.2002  0.0338]7  (4.103) 


The  objective  function  value  (4.59)  for  this  vector  is 

C(xopt)  = 0.4683 

The  Youla  parameter  corresponding  to  x°Pl  is 

-0.20022  + 0.0338 


Qi  (z)  := 


0.0008 


(4.104) 


(4.105) 


The  final  form  (equation  4.58)  of  the  robust  predictive  controller  for  this  Youla  parameter  is 


Q(z)  = 


0.0798^+0.03372  -0.0135 
22+ 0.00022 -0.1682 


(4.106) 


The  modified  design  procedure,  using  equations  (4.85)  and  (4.86)  in  the  optimization 
(4.55)  - (4.56),  yields  the  optimal  vector 


xopt=  [0.2448  -0.0959  -0.0332]1 


(4.107) 


which  corresponds  to  the  following  Youla  parameter 

— 1"\  —0.09592  - 0.0332 


Qli(z)  = 


\ z J 


2 + 0.2448 


(4.108) 


The  final  form  (equation  4.58)  of  the  robust  predictive  controller  for  this  Youla  parameter  is 


cu(z)  = 


0.184l23+0.089622  -0.01152-0-0133 
23  - 0.275722  - 0.55832  - 0. 1660 


(4.109) 


The  servo  performance  of  the  nominal  predictive  controller,  along  with  both  robust 


controllers,  is  shown  in  Figure  4.10.  No  uncertainty  is  assumed  present  in  the  plant 
model,  i.e.,  8p=0.  A unit  step  change  in  the  set  point  occurs  at  t = 0 and  t = 50  sec.  A 
step  disturbance  of  magnitude  0.2  occurs  at  t = 12  sec.  By  design,  both  the  nominal 
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Figure  4.10  Performance  comparison  of  the  three  different  control  designs:  the  nominal 
controlled:),  robust  controller  without  integral  action  (-.)  and  with  integral 
action 
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Figure  4.11  Control  actions  produced  by  the  three  different  control  designs:  the  nominal 
controlled:),  robust  controller  without  integral  action  (-.)  and  with  integral 
action  (--). 
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Figure  4.12  Performance  comparison  of  the  different  control  designs  for  the  perturbed 
plant:  the  nominal  controlled:),  robust  controller  without  integral  action  (-.) 
and  with  integral  action 


Figure  4.13  Control  actions  produced  by  the  different  control  designs  for  the  perturbed 
plant:  the  nominal  controlled:),  robust  controller  without  integral  action  (-.) 
and  with  integral  action  (— ). 
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predictive  controller  and  the  robust  controller  with  integral  action  exhibit  no  steady-state 
offset,  while  the  robust  controller  designed  without  integral  action  does  produce  offset. 
The  disturbance  rejection  properties  of  the  three  designs  are  quite  similar. 

Figure  4.1 1 shows  the  control  actions  produced  by  the  three  control  designs.  Again, 
all  three  controllers  show  similar  behavior,  and  it  is  noted  that  all  three  designs  require  only 
modest  control  actions  in  regulation  and  disturbance  rejection  for  the  nominal  plant. 

Figure  4.12  shows  the  response  of  the  three  predictive  controllers  when  the  plant  is 
given  by  the  model  (4.101).  As  noted  above,  the  nominal  controller  produces  an  unstable 
closed  loop  for  this  specific  plant  model,  and  the  vertical  range  has  been  restricted  to 
highlight  the  responses  of  the  two  stable  designs.  The  robust  controllers  both  produce 
stable  closed  loops  with  the  perturbed  plant,  however,  there  are  more  noticeable  transients 
after  the  set  point  changes.  The  disturbance  rejection  properties  of  the  two  robust  designs 
are  still  quite  good.  The  steady-state  offset  of  the  original  robust  control  design  is  still  very 
apparent. 

Figure  4.13  shows  the  control  actions  produced  by  the  three  control  designs  in 
controlling  the  perturbed  plant.  The  nominal  controller  requires  actions  of  larger  magnitude 
than  shown  on  this  plot;  the  vertical  range  is  truncated  to  highlight  the  control  actions  of  the 
robust  controllers. 

It  is  of  interest  to  find  out  if  it  is  possible  to  design  a predictive  controller  for  the 
nominal  plant  that  is  robustly  stable.  Table  4.1  shows  the  results  of  designing  predictive 
controllers  for  different  values  of  Ny  and  Nu.  The  value  of  % is  given  for  the  choices  of 
Ny  and  Nu  shown,  and  as  discussed  in  Chapter  2,  values  less  than  one  correspond  to  a 
robustly  stable  controller.  For  all  the  entries  of  Table  4.1,  the  move-suppresion  parameter 
A.  is  set  equal  to  zero.  For  all  the  entries  where  Nu  > 1 the  predictive  controller  given  in 
equations  (4.96)-(4.98)  is  recovered,  which  is  not  robustly  stable. 
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Table  4.1  Nyquist  robust  stability  margin  for  nominal  predictive  controllers 
designed  using  various  Ny  and  Nu  values. 


Ny  \ NU 

1 

2 

3 

4 

5 

2 

0.9165 

1.0954 

3 

0.8595 

1.0954 

1.0954 

4 

0.8332 

1.0954 

1.0954 

1.0954 

5 

0.8187 

1.0954 

1.0954 

1.0954 

1.0954 

4.7.2  Example  2 


The  second  example  is  the  stable  nominal  plant  model 

0.30z  + 0.42 
1 z1  + 1.20z  +0.54 


(4.110) 


with  the  corresponding  nominal  plant  parameter  vector  p°  =[1.20  0.54  0.30  0.42]'.  The 


matrix  Qp  for  this  example  is 


0.1000 

0.0490 

0.0140 

0.0105 

0.0490 

0.0700 

0.0275 

0.0130 

0.0140 

0.0275 

0.1200 

-0.0208 

0.0105 

0.0130 

-0.0208 

0.0600 

(4.111) 


A nominal  predictive  controller  is  designed  for  the  system  using  the  tuning  parameters 
Ny  = 4,  Nu  = 2,  and  X = 0,  resulting  in  the  following  controller  and  prefilter  polynomials 

R(z)  = z2  - 0.5402z  - 0.4598  (4.112) 

S(z)  = 0.2197z2  + 1.5860z  + 0.591 1 (4. 1 13) 

T(z)  = 1.3037z4  - 0.88 14z3  + 0.5 101z2  + 1 .4643z  (4. 1 14) 


The  nominal  closed  loop  polynomial  is 

A*(z)  = z4  + 0.7257z3 


(4.115) 
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and  the  polynomials  Ai(z)  and  A2(z)  are  taken  as  Aj(z)  = z2+0.7257z  and  A2(z)  = z2. 
Again,  the  pole  of  the  basis  function  is  taken  to  be  a = 0.  It  is  straightforward  to  show  that 
the  perturbation 

5p  = [0.07  -0.10  -0.27  0.10]T  (4.116) 

yields  8pTQ~'Sp  = 0.9627  and  is  therefore  allowable.  The  plant  model  resulting  from  this 
perturbation  is 


P(z)  = 


0.03Z  + 0.52 
z2  + 1.27z  + 0.44 


(4.117) 


The  closed  loop  polynomial  resulting  from  this  plant  model  and  the  nominal  controller  is 

G(z)  = z4  + 0.7364z3  - 0.5440z2  + 0.0208z  + 0. 105 1 (4. 1 18) 


which  has  a root  outside  the  unit  circle,  thus  the  nominal  controller  is  not  robustly  stable. 

Robust  controllers  with  and  without  integral  action  are  designed  for  the  choice  N = 2. 
The  optimal  vector  x for  the  robust  control  design  without  integral  action  is 

xopt  =[0.6397  0.3024  -0.4932  -2.6443  -1.7309]1  (4.119) 

The  objective  function  value  (4.57)  for  this  vector  is 


C(xopt)  = 0.8793 


The  Youla  parameter  corresponding  to  x°P'  is 


Q2(z) 


-0.4932z2-2.6443z- 1.7309 
z2  +0.6397z  + 0.3024 


(4.120) 


(4.121) 


The  final  form  (equation  4.56)  of  the  robust  predictive  controller  for  this  Youla  parameter  is 


C2(Z) 


-0.2735z5  - 1.3502z-2.2453z3  - 1.4338z2  -0.1335z  +0.1297 
z5  + 0.9732z4  + 0.5697z3  + 0.8074z2  + 0.2559z  - 0. 1009 


(4.122) 


The  modified  design  procedure,  using  equations  (4.72)  and  (4.73)  in  the  optimization 
(4.53)  - (4.54),  yields  the  optimal  vector 

xopt  =[0.5843  0.3764  0.0660  -1.0629  -0.8583]1  (4.123) 

which  corresponds  to  the  following  Youla  parameter 


85 


fz-l') 0.0660z2  - 1 .0629 z - 0.8583 


(4.124) 


{ z J z2  + 0.5843z  + 0.3764 


The  final  form  (equation  4.56)  of  the  robust  predictive  controller  for  this  Youla  parameter  is 


Note  that  the  degree  of  this  controller  is  only  5,  since  two  pole-zero  cancellations  at  the 
origin  occur  when  the  controller  is  constructed. 

The  servo  performance  of  the  nominal  predictive  controller,  along  with  both  robust 
controllers,  is  shown  in  Figure  4.14.  No  uncertainty  is  assumed  present  in  the  plant 
model,  i.e.,  8p=0.  A unit  step  change  in  the  set  point  occurs  at  t = 0 and  t = 50  sec.  A 
step  disturbance  of  magnitude  0.2  occurs  at  t = 12  sec.  By  design,  both  the  nominal 
predictive  controller  and  the  robust  controller  with  integral  action  exhibit  no  steady-state 
offset,  while  the  robust  controller  designed  without  integral  action  does  produce  offset. 
The  disturbance  rejection  properties  of  the  three  designs  are  quite  similar. 

Figure  4.15  shows  the  control  actions  produced  by  the  three  control  designs.  Again, 
all  three  controllers  show  similar  behavior,  and  it  is  noted  that  for  this  example  all  three 
designs  require  more  substantial  control  actions  in  regulation  and  disturbance  rejection  for 
the  nominal  plant. 

Figure  4.16  shows  the  response  of  the  three  predictive  controllers  when  the  plant  is 
given  by  the  model  (4.105).  The  nominal  controller  produces  an  unstable  closed  loop  for 
this  specific  plant  model,  so  the  vertical  range  has  been  restricted  to  highlight  the  responses 
of  the  two  stable  designs.  The  robust  controllers  both  produce  stable  closed  loops  with  the 
perturbed  plant,  however,  there  are  very  noticeable  oscillations  after  the  set  point  changes. 
The  steady-state  offset  of  the  original  robust  control  design  is  still  very  apparent. 


0.2857z5  + 0.8241z4  + 1.7302z3  + 2.5980z2  + 2.0468z  + 0.6249 
z5  + 0.7500z4  - 0.056  lz3  - 0.3488z2  - 0.8590z  - 0.4861 


(4.125) 


86 


1.4 


Q- 

D 

o 


1.2 

1 

0.8 

0.6 

0.4 

0.2 


IMA  . I 

I It  I \ A / 

II  II  u » 


0 


-0.2 


0 


20 


40  60 

Time,  t 


80 


100 


Figure  4.14  Performance  comparison  of  the  three  different  control  designs:  the  nominal 
controlled:),  robust  controller  without  integral  action  (-.)  and  with  integral 
action  (— ). 


Figure  4.15  Control  actions  produced  by  the  three  different  control  designs:  the  nominal 
controlled :),  robust  controller  without  integral  action  (-.)  and  with  integral 
action  (— ). 
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Figure  4.16  Performance  comparison  of  the  different  control  designs  for  the  perturbed 
plant:  the  nominal  controller(:),  robust  controller  without  integral  action  (-.) 
and  with  integral  action 


Figure  4.17  Control  actions  produced  by  the  different  control  designs  for  the  perturbed 
plant:  the  nominal  controller(:),  robust  controller  without  integral  action  (-.) 
and  with  integral  action  (— ). 
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Figure  4.17  shows  the  control  actions  produced  by  the  three  control  designs  in 
controlling  the  perturbed  plant.  The  nominal  controller  requires  actions  of  larger  magnitude 
than  shown  on  this  plot;  the  vertical  range  is  truncated  to  highlight  the  control  actions  of  the 
robust  controllers.  The  robust  design  with  integral  action  requires  smaller  initial  control 
actions  than  the  standard  robust  design.  Both  designs  require  similar  control  efforts  after 
the  onset  of  the  disturbance. 

Table  4.2  shows  the  values  of  the  Nyquist  robust  stability  margin  for  several  different 
predicitive  controllers  designed  using  the  values  of  Ny  and  Nu  shown.  For  all  the  designs, 
X = 0.  It  is  interesting  to  note  that  as  with  Example  1,  only  the  designs  with  Nu  = 1 are 
robustly  stable.  However,  the  Nyquist  robust  stability  margin  does  not  continually 
decrease  with  increasing  Nu,  as  is  seen  from  the  first  column  of  Table  4.2.  The  entries 
marked  with  an  X represent  Ny  and  Nu  pairs  that  result  in  a controller  that  is  not  nominally 
stable. 


Table  4.2  Nyquist  robust  stability  margin  for  nominal  predictive  controllers 
designed  using  various  Ny  and  Nu  values. 


Ny  ^ NU 

1 

2 

3 

4 

5 

2 

0.8767 

X 

3 

0.8372 

3.2163 

X 

4 

0.8701 

2.1186 

1.7620 

X 

5 

0.8690 

2.0701 

1.7132 

1.8914 

X 

4.8  Design  Equations  for  Nominal  Predictive  Control 

This  section  provides  specific  design  equations  used  to  synthesize  a nominal  predictive 
controller  following  the  approach  of  Crisalle  et  al.  (1989).  An  equivalent  formulation  is 
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given  by  McIntosh  et  al.  (1991).  The  final  design  equations  for  the  polynomials  (4.6)  - 
(4.8)  that  appear  in  the  predictive  control  law  (4.5)  are: 


R(z)  = zn 


,-i  Sr 


i+^'Xwr1) 

i'=l 


d-z"1) 


S(z) = zn 


'Ny  , 

r=l 


N, 


T (z)=lkiZl 


r=i 


(4.126) 


(4.127) 


(4.128) 


where  the  design  operators  F ,-(z_1)  and  T^z-1)  and  the  coefficients  kL,  i=  1,2,  ...,  Ny  are 
determined  from  the  process  model  according  the  following  procedure.  First,  rewrite  the 
nominal  plant  model  (4.1)  - (4.3)  in  the  equivalent  form 

A_j(z_1)  y(z)  = z_1  B_j(z-1)  u(z)  (4.129) 


involving  inverse  powers  of  z,  where  A_j(z_l)  and  B_j(z-1)  are  related  to  (4.2)  and  (4.3) 
in  an  obvious  manner  and  are  of  the  form 

A_i(z  1)  = l + r3_]jZ  z 2+...  + a_inz  n (4.130) 

B_!(z-1)  = b_lj0  +Z?_iiz-1  +...  + b_lnhz~nb  (4.131) 

To  obtain  the  design  operators  F^z-1),  which  are  polynomials  of  degree  n (the  order  of  the 
plant),  solve  the  set  of  Diophantine  equations 

Ei(z"1)A(z_1)A_1(z"1)  + z",  Fi(z_1)  = l,  i=  1,  2,  ....  Ny  (4.132) 

which  also  yields  the  intermediate  polynomials  E((z_1)  of  degree  M.  The  second  design 
operators,  the  polynomials  T,(z_1)  of  degree  are  obtained  by  decomposing  the 
product  Et(z-1)B_1(z-1)  in  the  form 

E,(z-1)B_1(z-1)  = G,(z-1)  + z-,r,(z_1) 


(4.133) 


90 


where  polynomials  G /(z-1)  of  degree  j'-l  are  known  as  the  dynamic  polynomials,  and  are 
characterized  by  the  fact  that  their  coefficients  are  the  sampled  values  of  the  step-response 
of  the  plant  (2.6.10).  Note  that  the  polynomials  T^z-1)  are  identically  zero  if  n\>  equals  0. 
In  turn,  the  coefficients  of  the  dynamic  polynomials  are  used  to  define  the  nonzero  elements 
of  the  Toeplitz  matrix  G^u  known  as  the  truncated  dynamic  matrix,  which  contains  only 
Nu  columns.  Finally,  the  coefficients  /q , i=\,  2,  ...,  Ny  are  obtained  as  the  components 
of  the  gain  vector  A:T=[k j k2  ...  k^y],  calculated  from  the  expression 


A design  technique  to  synthesize  robust  predictive  controllers  for  systems  subject  to 
ellipsoidal  uncertainty  has  been  presented.  The  technique  uses  a constrained  quasi-convex 
optimization  to  determine  the  coefficients  of  a fixed  order  Youla  parameter  to  robustly 
stabilize  the  system.  The  robust  controller  retains  the  nominal  servo  performance  of  the 
original  predictive  controller  designed  for  the  nominal  system.  A straightforward 
modification  is  given  to  allow  the  incorporation  of  integral  action  into  the  robust  controller 
design.  Furthermore,  the  design  can  accommodate  unstructured  uncertainty,  as  long  as  the 
problem  remains  equivalent  to  a rank  one  M-A  problem. 


(4.134) 


4.9  Conclusions 


CHAPTER  5 

ANALYTIC  SOLUTION  TO  A LIMITING  SYNTHESIS  PROBLEM 


5.1  Introduction 

This  chapter  presents  an  analytic  solution  to  the  limiting  case  of  the  robust  control 
design  problem  discussed  in  Chapter  2.  The  analysis  method  developed  in  Chapter  2 is  not 
strictly  applicable  to  the  limiting  case  of  a first  order  plant  and  constant  controller.  This  is 
the  case  where  the  characteristic  polynomial  is  of  degree  1,  i.e.,  k :=  n + tn  = 1.  For  this 
case,  an  analytic  solution  to  the  robust  synthesis  problem  is  available. 

This  chapter  is  organized  in  the  following  manner.  Section  2 presents  the  limiting  case 
discrete  time  problem  definition.  Section  3 presents  the  results  of  the  robustness  analysis 
for  the  special  case.  These  results  are  used  to  produce  the  solution  to  the  robust  synthesis 
problem,  which  is  presented  in  Section  4.  The  continuous  time  limiting  case  robustness 
design  is  presented  in  Section  5.  An  example  is  presented  in  Section  6 that  illustrates  the 
different  possible  results  of  the  robust  synthesis  procedure.  Concluding  remarks  are 
presented  in  Section  7. 


5.2  Problem  Statement 

Consider  the  closed  loop  system  shown  in  Figure  5.1,  where  the  compensator  is  a 
static-gain  controller  k 6 9?,  with  k * 0,  and  the  plant  is  the  discrete  time  system 

Pfcp  >-r^r  <s-‘> 

which  is  characterized  by  the  parameter  vector  p = [a  (3]T. 
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L 

Figure  5.1.  Feedback  loop  with  plant  P(z)  and  controller  k. 


The  plant  parameter  vector  is  modeled  as 

p = p°  + 5p  (5.2) 

where  p°  = [a0  (3°]T , (3°  * 0,  is  a vector  containing  the  known  nominal  values  of  the  plant 
parameters,  and  the  uncertainty  8p  is  an  element  of  the  set 

Ep  :=  {5p  e 9*2  | Sp^-'Sp  < l}  (5.3) 


9*2x2 


where 

q 9il  912 

L 912  922 j 

is  a symmetric  positive-definite  matrix. 

With  the  plant  model  (5.1)  placed  in  the  closed-loop  of  Figure  5.1,  the  characteristic 
polynomial  is  simply  G(z)  = z + g where  g = a + k(3.  It  is  obvious  that  the  closed  loop 
will  be  robustly  stable  if  and  only  if  all  possible  values  of  g satisfy  Igl  < 1.  Therefore,  it 
must  be  established  how  the  ellipsoidal  uncertainty  (5.3)  affects  the  characteristic 
polynomial  coefficient. 

The  relationship  between  the  plant  parameter  vector  p and  the  characteristic  polynomial 
coefficient  g is  given  by  the  linear  map  g = Sp  where  S = [1  k].  The  nominal  values  of 
the  plant  parameters  map  to  a nominal  characteristic  coefficient  g°  = Sp°  or,  explicitly, 
g°  = a0  + k(3°.  It  follows  that  the  uncertainty  in  the  coefficient  of  the  characteristic 
polynomial  8g  = g - g°  is  of  the  form 

Sg  = S8p  (5.4) 

where  Sg  is  a real  scalar.  As  8p  takes  on  all  values  in  Ep,  (5.4)  generates  a family  of 
characteristic  polynomials 

g = { G(z)  I g = g°  + Sg,  Sg  = SSp,  Sp  g Ep } 


(5.5) 
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It  is  not  difficult  to  show  that  when  8p  belongs  to  the  set  Ep  defined  in  (5.3),  then  5g 
given  by  (5.4)  lies  in  a continuous,  real  interval  Eg  that  is  symmetric  with  respect  to  zero. 
That  is,  Eg  is  of  the  form 

Eg  = {Sg  e 91  | |8g|  < Sg*,  Sg*e  91,  Sg*>  0}  (5.6) 

The  general  theory  of  Guzzella  et  al.  (1991)  is  not  applicable  to  the  limiting  plant/controller 
pair  considered  in  this  chapter  because  the  line-segment  (5.6)  is  in  fact  a degenerate 
ellipse.  Since  the  condition  for  robust  stability  is  Igl  < 1,  the  endpoint  8g*  of  the 
characteristic  interval  Eg  is  of  particular  significance  to  the  analysis  of  robustness.  The 
explicit  form  for  Sg*  is  given  in  Theorem  5.1 

Theorem  5.1  Consider  the  closed-loop  system  of  Figure  5.1,  with  k e 9?  and  the 
family  of  plants  P(z)  described  by  (5.1)-(5.3).  Then  the  uncertainty  8g  in  the  coefficient  of 
the  characteristic  polynomial  lies  in  the  interval 

£g  = j Sg  e 91  | |8g|  < ^|sQSf]|  (5.7) 

Proof.  Writing  the  endpoint  of  the  interval  Eg  in  (5.6)  as  8g*  = S8p*  where  8p* 
some  element  of  Ep,  it  follows  by  linearity  that  8p*  must  lie  on  the  boundary  of  the  ellipse 
Ef,  i.e.,  Sp*  satisfies 

8p*TQ“'Sp*  = 1 (5.8) 

The  value  of  8g*  is  found  by  solving  the  optimization  problem 

Sg*  = max  S8p  (5.9) 

8pe  £p 

where  Ep  is  defined  in  (5.3).  This  constrained  optimization  problem  is  solved  using  the 
Lagrange  multiplier  method.  Let 

J(8p)  = SSp  - ^(SpTQ-18p  - l) 

be  the  Lagrangian,  where  X > 0 is  the  unknown  multiplier.  Equating  to  zero  the  derivative 
of  J(Sp)  with  respect  to  8p  and  solving  for  the  optimal  value  of  Sp  yields 

Sp*  = r‘QST 


(5.10) 
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The  Lagrange  multiplier  is  readily  obtained  by  substituting  (5.10)  into  the  constraint  (5.8), 
yielding  \2=  SQST.  Substitution  of  (5.10)  into  5g*  = S8p*  yields 

8g*  = r'SQST  (5.11) 

from  which  it  follows  that 

5g*  = (SQST  )1/2  (5-12) 

The  quadratic  form  SQST  = qn  + 2qi2k  + q22^2  is  a positive  scalar  for  all  values  of  k, 
since  S is  never  equal  to  the  null  vector  and  Q is  a positive-definite  matrix.  Therefore,  the 
uncertainty  Sg  satisfies 

|8g|  < VSQST  (5.13) 

and  the  uncertainty  region  2^  given  in  (5.6)  can  be  rewritten  in  the  form 

2g  = { 5g  e 9t  I |8g|  < t/sQS7  } V 

It  is  of  interest  to  note  that  if  the  continuous-time  limiting  case  is  considered,  the 
resulting  closed  loop  characteristic  polynomial  is  exactly  identical  to  the  discrete  time  case 
except  that  the  discrete  variable  z is  replaced  by  the  Laplace  variable  s.  Therefore,  the 
analysis  of  the  previous  section  holds  without  modification  except  in  this  case  the  robust 
stability  criterion  is  that  g be  positive  for  all  values  of  the  uncertainty  5g.  The  continuous 
time  case  will  be  considered  further  in  Section  5. 

5.3  Robustness  Analysis 

This  section  presents  a criterion  for  testing  whether  a given  controller  k places  the 
roots  of  all  of  the  members  of  the  characteristic  polynomial  family  £ inside  the  unit  circle. 

Theorem  5.2  Consider  the  closed-loop  system  of  Figure  5.1,  with  k e 9^  and  the 
family  of  plants  P(z)  described  by  (5.1)-(5.3).  Furthermore,  let  q = SQST.  Then  the 
closed  loop  system  is  robustly  stable  if  and  only  if 

lg°l  + Vq  < 1 


(5.14) 
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Proof.  Again,  the  closed  loop  is  robustly  stable  if  and  only  if  Igl  < 1.  Since 
g = g°  + 8g,  if  lg°  + 5gl  < 1,  then  the  theorem  is  proved.  From  the  triangle  inequality 
and  (5.13) 


|g|  = |g°  + &g 
By  the  definition  of  q,  it  follows  that 


< 


g 


^[sQS1 


(5.15) 


|g|  * 


It  is  not  hard  to  verify  that  equality  will  hold  for  some  member  of  the  uncertainty  set  (5.3). 
Therefore,  the  condition  Igl  < 1 is  equivalent  to 


g 


Vq  < 1 


Remark  1.  It  is  noted  that  the  inequality  lg°l  + Vq  < 1 is  a more  restrictive  condition 
than  lg°l  < 1,  so  a controller  must  be  nominally  stabilizing  to  be  robust. 

Remark  2.  The  condition  lg°l  + Vq  < 1 is  equivalent  to  checking  the  Schur  stability  of 
the  extreme  polynomials  of  g,  namely  Gj(z)  = z + g°  + Vq  and  G2(z)  = z + g°  - Vq  . This 
is  the  method  suggested  in  Barmish  (1994)  for  analyzing  the  robust  stability  of  low-order 
interval  polynomials. 

Theorem  5.2  provides  a necessary  and  sufficient  condition  for  testing  the  robustness 
of  a controller  k.  However,  a slightly  modified  version  of  this  condition  is  used  in  the 
synthesis  discussion  to  follow.  Therefore,  before  the  robust  synthesis  problem  is 
addressed,  a robustness  parameter  r(k)  is  defined  and  the  condition  presented  in  Theorem 
5.2  is  restated  in  terms  of  this  new  parameter.  Let 

r(k)  :=  1 - lg°l  - Vq  (5.16) 

where  r(k)  G 91  can  be  interpreted  as  the  difference  between  the  magnitude  of  the  maximal 
root  of  Q and  the  critical  magnitude  Igl  = 1 . 

Corollary  5.1.  The  feedback  loop  of  Figure  5. 1 is  robustly  stable  if  and  only  if  r(k)  > 0. 
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5.4  Robust  Synthesis 


This  section  considers  the  problem  of  finding  a controller  k that  yields  a robustly  stable 
closed  loop.  A robust  controller  must  yield  a nominal  characteristic  polynomial  that  is 
stable.  Let  7C  = {k  e 91  : lg°l  < 1}  be  the  set  of  all  controllers  that  lead  to  a stable 
nominal  closed  loop.  It  is  straightforward  to  show 

ac=  jk  e 

and 

X=  jk  € 

Using  the  result  of  Corollary  5.1,  define  the  (possibly  empty)  set  of  robustly 
stabilizing  controllers  Tfo  = (k  e 91  : r(k)  >0).  The  goal  of  the  robust  synthesis 
procedure  is  to  find  the  optimal  controller 

k = arg  max  r(k)  (5.18) 

ke9C 


-(1  + a ) 1-a 


P 


0 ’ o0 


PL 


if  P°  > 0 


(5.17a) 


^l-a°  -(1  + aV 

Q0  ’ r.0 


if  (3°  < 0 


(5.17b) 


The  controller  k is  called  optimally  robust  if  r(k)  >0.  To  highlight  the  dependence  on  the 
controller  parameter  k,  the  function  r(k)  can  be  written  in  the  form 

r(k)  = 1 - a0  + k(3°  - ^qn  + 2q,2k  + q22k2  (5.19) 


Although  (5.19)  is  a simple  function,  due  to  the  presence  of  the  absolute-value  operator  the 

derivative  with  respect  to  k does  not  exist  at  the  point 

ryO 

ko=-j^  (5.20) 

Consequently,  the  maximum  of  r(k)  cannot  be  computed  simply  by  setting  the  first 
derivative  of  r(k)  with  respect  to  k equal  to  zero  and  solving  for  an  optimal  k.  This 
problem  is  avoided  using  an  alternate  definition  of  r(k).  Let 

r(k)  = min  { n(k) , r2(k)  } (5.21) 

where 


and 


r i (k)  :=  1 - g°  - Vq 


(5.22) 
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r2(k)  :=  1 + g°  - Vq 


(5.23) 


It  follows  that  r(k)  = ri(k)  when  g°  > 0,  and  r(k)  = r2(k)  when  g°  < 0.  Furthermore, 
ri(ko)  = r2(ko)  since  g°(ko)  = 0,  where  ko  is  defined  in  (5.20).  Note  that  rj(k)  and  r2(k) 
are  continuously  differentiable  functions  of  k. 

The  problem  of  finding  the  maximum  of  r(k)  is  solved  by  analyzing  individually  the 
maxima  of  rj(k)  and  r2(k).  An  analytical  solution  is  possible  because  functions  rj(k)  and 


r2(k)  are  concave.  Denote  as 


k,  = - 


(3° 


91.922  9l2 


0x2 


9 22  - (P°) 


dr  i 


the  solution  to  ^ = 0,  and  denote  as 


k2  = 

922 


922  922 

P |9ll922  ~ 9 1 2 


(5.24) 


0x2 


922  \ 922  (P  ) 


(5.25) 


dr? 

the  solution  to  ^7  = 0.  If  k]  and  k2  are  real,  they  are  the  global  maximizers  of  ri(k)  and 


r2(k),  respectively.  The  numerator  of  the  fraction  inside  the  radical  in  (5.24)  and  (5.25)  is 
positive  because  q1i922-q12  = det(Q)  > 0.  If  q22-(P0)2  < 0,  then  both  k\  and  k2  are 
complex,  and  neither  of  the  branch  functions  ri(k)  and  r2(k)  has  a real  maximum.  Theorem 


5.3  presents  a comprehensive  solution  to  the  robust  synthesis  problem  in  terms  of  the  three 
controllers  ko,  k),  and  k2. 


Theorem  5.3.  Consider  the  closed-loop  system  of  Figure  5.1,  with  k e 91  and  the 
family  of  plants  P(z)  described  by  (5.1)-(5.3).  Also  define  ko,  k)  and  k2  as  in  (5.20), 


(5.24),  and  (5.25),  respectively.  Then  the  optimal  controller  is  given  by 


k =i 


ki 

k2 


l k0 


if  a0  + k 1 [5 0 > 0 and  kj  e 91 
if  a0  + k2 P 0 < 0 and  k2  e 91 
otherwise 


Furthermore,  if  r(k)  > 0,  then  the  optimal  controller  is  robust. 


(5.26a) 

(5.26b) 

(5.26c) 


Proof.  The  proof  of  Theorem  5.3  is  given  in  the  Appendix. 
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Theorem  5.3  presents  the  complete  analytic  solution  to  the  limiting  case  robust 
synthesis  problem  of  a first  order  plant  subject  to  ellipsoidal  uncertainty  and  constant 
controller.  The  optimal  controller  is  found  simply  by  checking  the  sign  of  g°  at  the 
controllers  kj  and  k2-  If  neither  of  these  two  controllers  satisfies  the  conditions  (5.26a) 
and  (5.26b),  then  the  optimal  controller  is  ko-  It  is  remarked  that  the  algorithm  in  Theorem 
5.3  results  in  a unique  choice  of  controller. 

5.5  Continuous  Time  Limiting  Case 

In  the  limiting  case  for  continuous  systems,  the  plant  (5. 1)  is 

Kr-v'-jh  (5'27) 

The  uncertainty  in  the  coefficients  is  still  given  by  (5.2)-(5.3),  and  the  controller  is  still  a 
static  gain  k.  The  resulting  characteristic  polynomial  is  G(s)  = s + g where  g = g°  + 8g, 
and  the  uncertainty  5g  lies  in  the  interval 

£g  = { 5g  e 51  | |8g|  < VSQST  } (5.28) 

For  the  continuous  time  case,  the  necessary  and  sufficient  condition  for  robust  stability  is 
that  the  characteristic  polynomial  coefficient  g be  positive.  Since 

|Sg|  < VSQST  (5.29) 

it  follows  that  if  g°  is  less  than  or  equal  to  zero,  the  closed  loop  is  not  robustly  stable  since 
the  allowable  perturbation  8g  = 0 yields  an  unstable  closed  loop.  Therefore,  it  is  assumed 
that  g°  > 0,  so  that  the  nominal  closed  loop  is  stable.  In  this  case,  the  only  way  that  the 
coefficient  g can  become  negative  is  by  passing  through  the  origin,  implying  that  g = 0 for 
an  allowable  8g.  This  implies  that  if  the  nominal  system  is  stable,  then  robust  stability  is 
guaranteed  if  and  only  if  g * 0 for  all  allowable  8g.  This  condition  is  equivalent  to 

g°  > VSQST  (5.30) 

For  the  continuous  time  case,  the  robustness  parameter  r(k)  is  defined  as 
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r(k)  = g°  - Jioi7 


(5.31) 


Corollary  5.2.  The  continuous  time  limiting  case  feedback  loop  is  robustly  stable  if  and 
only  if  r(k)  > 0. 

Computing  the  range  of  values  of  k that  lead  to  a stable  nominal  closed  loop  is 
straightforward.  Let  3C={ke  SK  : g°  > 0 } be  the  set  of  all  controllers  that  lead  to  a 
stable  nominal  closed  loop.  It  is  straightforward  to  show 


Using  the  result  of  Corollary  5.2,  define  the  possibly  empty  set  of  robustly  stabilizing 
controllers  = {k  e !?Cl  r(k)  >0}.  The  goal  of  the  robust  synthesis  procedure  is  to 
find  the  optimal  controller 


This  controller  is  optimal  only  if  two  conditions  are  met:  1)  a0  + ki(3°  > 0 and  2)  k\  is  real. 

If  (P0)2  > q22  then  ki  is  imaginary,  and  furthermore  r(k)  — > «=  as  |k|  — » °°  {where 
the  sign  of  k is  chosen  so  that  g°  > 0}.  In  this  case  there  is  no  finite  controller  that 
maximizes  the  robustness  parameter  r(k),  and  an  appropriate  robustly  stabilizing  controller 
should  be  chosen  using  additional  criteria.  This  situation  is  not  as  simple  as  the  discrete 
time  case  where  the  optimum  controller  can  be  found  directly.  However,  in  this  case  a 


(5.32a) 


and 


(5.32b) 


kopt  = arg  max  r(k) 


(5.33) 


Assuming  that  g°  > 0 then 


(5.34) 


(5.35) 
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robustly  stabilizing  controller  is  guaranteed  to  exist,  even  if  a large  (magnitude)  value  of  k 
is  required  to  stabilize  the  loop. 

If  (P0)2  < q22  then  k[  is  real.  This  implies  that  if  k\  satisfies  a0  + kiP°  > 0,  then  it  is 

the  unique  solution  to  (5.33),  i.e.,  ki  = kopt.  If  a0  + ki  (3°  < 0,  then  there  is  no  maximum. 

5.6  Example 

The  theory  discussed  in  this  chapter  is  illustrated  with  a simple  example.  Consider  the 
nominal  plant  P(z;  p°)  = where  P°  = 0.08  and  a0  will  take  different  values  for  this 
example.  Let  the  parametric  uncertainty  description  £^be  described  by  (5.3)  through  the 
matrix 

[ 0.02  -0.01  ‘ 

Q = L -0.01  0.02  . 

When  a0  = 0.2,  condition  (5.26a)  of  Theorem  5.3  is  satisfied  because  kj=-0.094e  01 
and  a°+k]P°  = 0.192  > 0;  hence,  the  global  maximizer  of  r(k)  is  k = -0.094.  The 
branches  rj(k)  and  r2(k)  are  plotted  along  with  the  function  r(k)  in  Figure  5.2a.  Note  that 
in  this  case,  both  branch  maxima  kj  and  k2  lie  to  the  right  of  the  intersection  point  k0. 
When  a0  = -0.25,  condition  (5.26b)  of  Theorem  5.3  is  satisfied  because  k2=1.094e  01 
and  a°+k2P°  = -0.162  < 0;  hence,  the  global  maximizer  of  r(k)  is  k = 1.094.  The 
functions  r(k),  rj(k),  and  r2(k)  are  plotted  in  Figure  5.2b,  where  in  this  case  both  branch 
maxima  lie  to  the  right  of  k0.  When  a0  = -0.04,  condition  (5.26c)  of  Theorem  5.3  is 
satisfied  because  o^+kjP0  = -0.048  < 0 and  a°+k2P°  = 0.048  > 0;  hence,  the  global 
maximizer  of  r(k)  is  k = Icq  = 0.50.  The  results  of  this  case  are  presented  in  Figure  5.2c, 
where  the  maximum  occurs  at  the  intersection  point  k = Icq. 
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(a)  (b) 


(c) 

Figure  5.2.  Plots  of  the  branch  functions  r^k),  r2(k),  and  r(k)  for  the  three  cases  of  the 
example.  The  global  maximizers  are:  (a)  k = k]  = -0.094,  (b)  k = k2  = 
1.094,  and  (c)  k = k0  = 0.50. 


5.7  Conclusions 

The  limiting  case  design  problem  considered  in  this  chapter  is  indeed  short,  simple, 
and  straightforward.  The  payoff  is  that  the  result  for  testing  robustness  is  simply  checking 
the  positivity  of  a function  containing  the  nominal  characteristic  coefficient  and  the  elements 
of  the  matrix  describing  the  uncertainty  in  the  plant  model  parameters.  Also,  the  procedure 
for  finding  the  maximally  robust  controller  (for  the  discrete  time  case)  is  a short  if-then-else 
sequence,  readily  adaptable  to  computer  implementation. 


The  case  where  maximization  leads  to  the  controller  kj  or  k2  (for  the  discrete  time  case 
- k]  for  the  continuous  time  case)  in  fact  appears  to  be  rather  rare.  To  see  why  this  is  so, 
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notice  that  to  satisfy  the  condition  leading  to  this  case,  the  variance  of  the  estimate  of  the 
process  gain  must  be  greater  than  the  square  of  the  estimate  itself.  The  examples  presented 
consider  a plant  with  a small  gain,  where  the  variance  of  the  gain  estimate  is  a large 
percentage  of  the  nominal  value. 

Although  the  above  discussion  appears  to  imply  that  the  results  presented  are  of  limited 
usefulness,  this  is  not  true.  It  so  happens  that  the  robust  stability  criteria  that  are  derived  in 
this  paper  are  the  same  criteria  that  occur  at  the  singular  frequencies  0 and  n in  the  general 
discrete  time  case.  The  scrutiny  provided  in  this  chapter  thus  provides  insight  into  the  more 
general  case. 


CHAPTER  6 
FUTURE  DIRECTIONS 

6.1  Ellipsoidal  Systems  with  Delays 

Many  systems  of  practical  interest  are  subject  to  delay.  In  the  continuous  time  case, 
the  presence  of  a delay  term  transforms  the  characteristic  polynomial  into  a so-called 
quasipolynomial.  The  results  of  Chapter  2 are  not  directly  applicable  to  this  delay  case.  In 
the  discrete  time  case,  if  the  delay  is  an  integer  multiple  of  the  sampling  time,  then  the 
resulting  system  is  still  a ratio  of  polynomials  and  the  results  presented  in  this  work  are 
applicable.  It  is  of  interest  to  see  if  stability  analysis  results  for  delay  ellipsoidal  systems 
can  be  derived.  As  a first  step,  a constant  delay  should  be  considered,  and  if  analysis 
results  are  obtained,  the  case  of  variable  delay  could  be  analyzed. 

6.2  Performance  Constraints 

Most  controllers  are  not  designed  solely  on  the  basis  of  stability  robustness. 
Performance  requirements  are  usually  present,  although  it  may  not  be  possible  to  directly 
incorporate  these  requirements  into  the  controller  design.  Usually  these  constraints  are 
posed  as  time-domain  constraints,  such  as  settling  time  limits,  but  they  may  also  be 
constraints  such  as  input  saturation.  If  a convex  optimization  can  be  constructed  to 
maximize  the  performance  of  the  system,  then  the  stability  constraint  that  is  optimized  in  the 
controller  design  methodology  of  Chapter  4 can  be  taken  simply  as  a constraint  in  this  new 
convex  optimization.  The  controller  design  procedure  then  becomes  a robust  performance 
problem. 
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6.3  LMI-based  Robust  Control  Design 

It  has  been  proposed  (Rantzer  1993)  that  the  finite  dimensional  optimization  for 
constructing  robust  controllers  can  be  reformulated  as  an  LMI  problem,  or  more 
specifically,  as  a generalized  eigenvalue  problem  subject  to  an  LMI  constraint.  Such  a 
formulation  would  be  able  to  take  advantage  of  the  powerful  methods  for  solving  LMIs  that 
have  been  developed  recently.  The  matrices  involved  in  the  LMI  are  constructed  from  a 
state  space  realization  of  the  transfer  function  vectors  T\  and  Tj  along  with  the  basis 
functions  tXk  and  [3k-  If  it  is  possible  to  derive  a method  for  constructing  these  state  space 
matrices  given  the  nominal  plant,  controller,  and  basis  functions,  then  this  optimization 
method  would  be  far  more  efficient  than  employing  standard  quasi-convex  optimization 
routines. 


APPENDIX 


1.  Proof  of  Lemma  2.2 


The  case  of  co  = 0 is  proved  first.  From  Figure  2.3  it  follows  that  pc(0)  is  half  of  the 
length  of  the  line  segment  Eo-  From  (2.23)  it  follows  that 

P c(°)  = V^7 

which  is  equation  (2.25b). 

For  to  > 0,  the  function  pc(to)  corresponds  to  the  distance  from  the  center  of  the  ellipse 
Eco  to  its  boundary,  and  can  therefore  be  considered  the  length  of  a vector  lying  on  the 
boundary  of  E^.  Let  this  vector  be  represented  by  St*((D)  which  satisfies 

(8t  * (co))T  Q”1  (8t  * (co))  = 1 (A  1.1) 

and  by  definition 

||St  * (o))||2  = pc(co) 

Considering  equation  (2.19),  Sx*(co)  also  represents  a perturbation  polynomial  AG(/C0*) 
that  can  be  written  in  Cartesian  and  polar  forms  as 

AGO'co*)  = 8tr(o))  + ;8xi(co)  = pc(o))e70  (A  1.2) 

for  some  angle  0.  Using  equations  (2.18)  and  (2.21),  the  positive  definite,  symmetric 
matrix  Qu  can  be  written  as 


Qco  = 


’A(to)  B(co)" 

_B(co)  D(to)_ 

wJ(co)Qgw^(co) 

wJ(a))Qgw^(o)) 


wJ(co)QgW/(co) 

wJ(C0)QgW7(G)) 


(A1.3) 


It  is  noted  that  the  elements  A(oo),  B(co),  and  D(co)  are  real  polynomials  in  to.  Equation 
(A  1.3)  and  the  definition  of  8x  can  be  used  to  rewrite  (Al.l)  as 
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D(8t^)2  + A(dx*j)2  - 2B(8t^)(8t/)  = det(Qw)  (A  1.4) 

where  the  dependence  on  co  has  been  dropped  from  all  terms.  Now,  from  (A  1.2)  it 
follows  that 

S'tfl  = pc(O))cos(0) 
and 

8i/  = pc(co)sin(0) 

Substituting  these  expressions  into  (A  1.4)  yields 

p2  (co)[d  cos2  (0)  + A sin2  (0)  - 2B  sin(0)cos(0)]  = det(Qu ) 
or 

Pc(®)  = 


*\2 


Vd  et(Qm) 


Vdcos2(0)  + Asin2(0)  - 2Bsin(0)cos(0) 


(A1.5) 


All  that  remains  is  to  determine  0.  Write  G°(/'co)  in  polar  form  as 


G(;‘co)  = rze 


where  from  (2.19)  it  follows  that 


rt  = 


(A  1.6) 


(A1.7) 


GuO’(o)  = ||t(cd)||2 

The  vectors  given  by  (A  1.2)  and  (A  1.6)  are  colinear,  however,  they  point  in  opposite 
directions.  This  implies  cos(a)  = -cos(0)  and  sin(a)  = -sin(0).  The  real  and  imaginary 
parts  of  G°(j(Q)  can  be  expressed  in  terms  of  the  quantities  in  equation  (A  1.6)  as 
xR  = rx  cos(a)  and  T/  = rTsin(a).  These  relationships  can  be  combined  to  give 


and 


cos(0)  = 


sin(0)  = 

rx 


Substituting  these  expressions  into  (A  1.5)  and  simplifying  yields 

rx  ydetlQJ 


Pc(“)  = 


d(t|)  + a(t5)-2B(twt/) 


Finally,  using  equations  (A  1.3)  and  (A  1.7)  results  in 
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Pr(CO)  = I 

^XT(03)Q“1X((0) 

which  is  equation  (2.25a).  This  completes  the  proof.  V 

2.  Proof  of  Lemma  3.1 

The  proof  of  this  Lemma  is  broken  down  into  three  separate  parts,  each  covering  a 
portion  of  the  claims  in  the  Lemma. 

Property  1.  h{ co)  is  a real-rational  function  of  co2,  and  thus  an  even  function  of  0). 
Furthermore,  the  numerator  degree  of  this  function  is  4/c-8  and  the  denominator  degree  is 
4*-4. 

Proof:  First,  it  is  shown  that  all  elements  of  the  right-hand  side  of  (3.5)  are  polynomials  in 
CO.  Consider  the  elements  of  the  matrix  Qw  given  in  (3.3).  From  the  form  of  the  vector  wr 
it  is  not  hard  to  see  that  A is  a polynomial  in  even  powers  of  co  with  a constant  term,  and 
the  degree  of  A is  2k-4  if  k is  even  and  2k-2  if  k is  odd.  Similarly,  D is  also  a polynomial 
in  even  powers  of  co,  but  with  no  constant  term.  The  degree  of  D is  2k-2  if  k is  even  and 
2k-4  if  k is  odd.  B is  a polynomial  in  odd  powers  of  co,  and  it  has  no  constant  term  either. 
The  degree  of  B is  2k-3.  Since  the  elements  of  Qg  are  real,  it  follows  that  all  these 
polynomials  have  real  coefficients.  Thus  the  numerator  of  h( co)  is  a polynomial  in  even 
powers  of  co  with  lowest  term  co2.  The  leading  degree  of  this  polynomial  is  2k-6.  The 
leading  coefficient  is  a minor  of  the  positive  definite  matrix  Qg  so  it  is  always  positive. 

Looking  at  the  definitions  of  Xr  and  X/ , it  is  straightforward  to  see  that  Xr  is  a 
polynomial  in  even  powers  of  co  with  a non-zero  constant  term.  The  degree  of  this 
polynomial  is  k if  k is  even,  and  k- 1 if  k is  odd.  Also,  it  can  be  seen  that  X/is  a polynomial 
in  odd  powers  of  co  with  no  constant  term,  with  degree  k- 1 if  k is  even,  and  k if  k is  odd. 
Thus  X£  and  Xj  are  polynomials  in  even  powers  of  co,  but  x^  has  a constant  term  while  x{ 
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does  not.  Also,  x^X/  is  a polynomial  in  odd  powers  of  co  with  no  constant  term.  From 
this  it  follows  that  all  three  elements  of  the  denominator  of  h{ co)  are  polynomials  (with  real 
coefficients)  in  even  powers  of  0)  with  lowest  term  to2.  It  is  straightforward  to  show  that 
the  degree  of  this  polynomial  is  Ak-2.  However,  both  the  numerator  and  denominator  have 
lowest  order  terms  to2.  Therefore,  a factor  of  co2  can  be  cancelled  from  each,  lowering  the 
numerator  degree  to  Ak-8  and  the  denominator  degree  to  Ak-A.  ♦ 

Property  2.  /z(c o)  is  a finite,  non-negative  function  for  co  > 0. 

Proof:  The  matrix  Qu  is  positive  definite  for  co  > 0.  Furthermore,  X is  never  equal  to  the 
null  vector.  Therefore,  xTQw*x  > 0 for  co  > 0.  From  (3.4)  it  follows  that  h~x{ co)  = 
xTQu'x.  Therefore,  h{ co)  < °°  for  co  > 0.  All  of  the  elements  of  qJ  are  bounded  for  finite 
co,  so  xTQw'x  < °°  for  0 < co  < °°.  Note  that  /z(co)  ->  0 as  co  ->  °°,  which  follows 
immediately  from  the  degree  conditions  of  Property  1 . Therefore, 

0 < h( co)  < « for  co  > 0 ♦ 


Property  3.  h( 0)  is  a finite,  positive  quantity. 

Proof:  In  the  proof  of  Property  1,  it  is  noted  that  both  numerator  and  denominator 
polynomials  of  h( co)  have  lowest  order  terms  co2.  Thus  h( 0)  = 0/0.  However,  once  the 
factor  of  co2  is  cancelled  from  both  polynomials  the  value  of  h( 0)  can  be  extracted.  It  is  a 
straightforward  but  tedious  exercise  to  show  that  this  process  yields 


H 0)  = 


4k,k(go)2 


2 

4k.k4k-l.k-l  ~ 4k,k-l 

+ 4k-i ,k-i (S? )2  ' 24k.k-i(go)(g?) 


(A2.1) 


where  q* j is  the  (i,j)  element  of  the  matrix  Qg.  Since  the  numerator  is  a minor  of  a positive 
definite  matrix,  it  is  a strictly  positive,  finite  quantity.  Thus 


V4k,k4k-l,k-l  > 4k,k-l 


(A2.2) 
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It  remains  to  show  that  the  denominator  of  (A2.1)  is  also  strictly  positive  and  finite.  The 
elements  qk,k  and  qk- 1 ,k-i  are  diagonal  elements  of  a positive  definite  matrix,  so  they  are 
positive.  Since  g®  and  gg  are  real  numbers,  it  follows  that  the  first  two  terms  in  the 
denominator  are  positive.  It  is  assumed  that  both  g°  and  gg  have  the  same  sign,  as  is 
necessary  if  the  nominal  closed  loop  is  to  be  stable.  Thus  if  q^k-i  ^ 0,  the  denominator  of 
(A2.1)  is  positive.  Now  consider  the  case  where  qk,k-i  > 0-  Using  (A2.2),  it  follows  that 

qk,k(gS)2  + qk-i,k-i(g?)2  - 2qk,k-i(go  g?) 


qk,k(go)2  + qk-l.k-ltg?)2  - 2Vqk,kqk-l,k-l  (go  g^) 


0x2 


0 „0x 


^k.k  (go)  ~ Vqk-l,k-l  (g?) 


> 0 


Since  this  proves  the  denominator  is  also  strictly  positive  and  finite,  the  final  result  is 

0 < h( 0)  < ♦ 


Taken  together,  the  three  properties  prove  all  the  claims  of  Lemma  3.1. 


3.  Proof  of  Theorem  3.2 


In  order  to  prove  the  Theorem,  the  relationships  between  A(co),  B(co),  D(co), 
and  T/(co)  and  A(j(o),  B(jco) , D(jco),  T/^jco),  and  ^/(j®)  are  derived  and  then  it  is 
shown  that  when  equations  (3.1 1)  and  (3.12)  are  substituted  into  equation  (3.10),  equation 
(3.5)  is  recovered.  The  first  step  is  to  write  the  various  polynomials  from  (3.5)  as 
polynomials  in  jco.  In  order  to  do  so,  the  vectors  wR  and  w /,  introduced  in  equation 
(2.18),  are  rewritten  as 


and 


Wp  = 


0 (j(0)4  0 (jco)4  0 (j(o) 


. xO 


T 


W, 


= or1  - (j®)5  o (jo,)4  o (jco)1  o 


• \3 


(A3.1) 
(A3. 2) 


It  follows  from  the  definition  of  T(co)  that 
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if  k is  even,  and 
if  k is  odd.  Defining 


**(®)  = g8+."  + g*-2(j®)*  2+(j“)A 

M®)  = go  + -.  + g*-i(j  “)*_1 


nic 


**(j®)  = X§2^(j«)2/ 

1=0 


where 


mR  = 


-k  k even 
2 

.^(fc-1)  ko  dd 


and  gk  :=  1 , it  follows  that 

M®)  = *r{  j®) 

Equations  (A3.3)-(A3.4)  are  simply  (3.1  la).  Similar  to  the  above  case, 

*/(®)  = 0'1)(g?(j®)1+...  + gJ-i(jto)*_1) 

if  k is  even,  and 

= Ci‘1)(gf(j«)1  +-+g°-i(jto)*"1)  + (j^Xj®)* 

if  k is  odd.  Defining 

m. 

0 


where 


and  gk  :=  1,  it  follows  that 


*/(j®)  = £g2*-i(j®) 

(=0 


\k  /ceven 

m,  = < , 2 

\Uk  + 1)  k odd 


*/(®)  = (T^/G®) 


or 


*/0®)  = jx/(a>) 


(A3. 3) 


(A3. 4) 


(A3. 5) 


(A3. 6) 


(A3 .7) 


(A3. 8) 


Equations  (A3.6)-(A3.7)  are  just  (3.1  lb). 


Ill 


A proof  of  all  equations  (3. 1 2)-(3. 1 3)  is  not  given,  rather,  one  polynomial  is 
discussed.  The  polynomial  A(co)  is  defined  in  equation  (2.34)  as 

A(g>)  = wJ(co)Qgw^((0) 

We  assume  k is  odd,  the  case  where  k is  even  is  similar.  Then 


wj (a>)Qgw*(a>)  = [(jco)*  1 •••  (jco)2  0 (jco)0]^,  •••  qk] 


k- 1 


(jco) 


(jto) 

0 

. (j©)°  . 


where  qj  represents  column  i of  the  matrix  Qg.  Thus 

A(co)  = [(jco)*'1  •••  Gco)°  ][qyt  Gco)°  + q^_2(ja))2  +•••  + Qi  Geo)*'1 
= (jco)°{qw(jco)0  +qw_2(jco)2  +...  + q/U(jco)*  '}  + •••  + 

Cjco)*-1  {qu  (jco)0  + %k-2  O®)2  + • • • + qu  Geo)*-1 } 

= Qk,k  +{(ik,k-2  )(jCO)2  + ...+ 

(q3,l  -+-qi,3)Gco)2A:_4  +qUiaco)2*-2 


Several  lengthy  calculations  allow  this  to  be  written  as 


1*  A 

A(co)  = A(jco)  = £a2^(jco)2£  (A3. 9) 

where  nA  = k- 1,  and  if  mA  is  defined  as  nA/2,  then 


i 

X 9(k-2h),(k-2C+2h)  0 < £ < mA 


a2r  - * 


h=0 


I 


9(k-2h),(k-2M-2h) 


(mA  + l)<^<nA 


(A3. 10) 


h=(C-mA) 

(A3. 9)  - (A3. 10)  are  simply  equations  (3.12a)  and  (3.13a).  For  the  case  where  k is  even, 
the  first  element  of  wr  is  zero,  so  the  highest  power  of  jco  is  2k-4  and  nA  = k- 2. 
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However,  the  formula  (A3. 10)  still  holds.  The  proofs  of  equations  (5.12b,c)  and 
(5. 13b,c)  are  entirely  similar  and  are  not  included.  However,  it  is  noted  that 

B(jco)  = jB(co)  (A3. 11) 

and 

D(jco)  = -D(to)  (A3. 12) 

where  B(co)  and  D(oo)  are  given  in  (2.34)  and  B(jco)  and  D(jco)  are  given  by  (3.12b)  and 
(3.12c).  Then,  substitution  of  (A3. 5),  (A3. 8),  (A3. 9),  (A3. 11),  and  (A3. 12)  into  (3.10) 

yields  equation  (3.5)  V 


4.  Proof  of  Theorem  3.3 


In  order  to  prove  the  Theorem,  the  relationships  between  A(co),  B(co),  D(co),  T/?(co), 
and  T/(co)  and  A(e7<u),  B(e7C0),  D(e7a)),  x^^03  j,  and  are  derived  and  then  it  is 

shown  that  when  equations  (3.27)  and  (3.28)  are  substituted  into  equation  (3.26),  equation 
(3.23)  results.  The  first  step  in  the  proof  is  expressing  the  vectors  xvR  and  w /,  as 


T 

xv  R = 


and 


[1(« 


j(k- l)ft)  + e~j(k- l)w 


) - ip +*-'") 


T 

W,  = 

' 1 U(k- 1)(0  -j(k- l)to\  . 

3 

•*-» 

1 

1 

3 

"**» 

0 

/ 

.2)1  1 

2 j\  1 

. 

(A4.1) 


(A4.2) 


From  the  definition  of  T(co),  it  follows  that 

Tfl(CD)  = wJ(co)g°  + cos(/cco) 


= ^eM-\)^+e-j(k- i)coj  ...  i 

= 7Ig?P  + e-*“) 


g°  + ^(ejk(0+e-jk'°) 


if  gk  :=  1.  Defining 


= 2x^(co) 


(A4.3) 
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then 


*,(«*)  = ls°l(ei'm  + <T*“) 


t=0 


which  is  equation  (3.27a).  In  similar  fashion 

t7(co)  = wJ(co)g°  + sin(/cco) 


i-p 

2;l 


7(/:-l)co  _ -j(k-\)m 


-e 


if  ^ :=  1 . Defining 
then 


) - ° 
= fig  ?p"  - 

Zj  e=o 

= 2/t,(0» 


g°  + 


(A4.4) 


= tg?(^-  - «-j“) 

^=0 


which  is  equation  (3.27b). 

A true  proof  of  all  equations  (3.28)-(3.31)  is  not  given,  but  the  formulation  of  the 
equations  for  the  polynomial  A(e;“)  in  Lemma  5.5  is  discussed.  The  polynomial  A(co)  is 
defined  in  equation  (2.34)  as 


A(co)  = wi(co)Q  Wfl(co) 


from  which  it  follows  that 


k k 

A(co)  = £ X(w*Mw/?)p(QA/ 

t= l p=i 


* k 


- 1 


j(k-p)  co  + e-j(k-p)  co 


*=I  p=l 
* * 


K, 


= i£  + ej(('p)a  + e~j(('p)(0  + e~j(2k'e'p)(0]qe  p 

i=\  p=  l 

It  is  noted  that  the  value  in  parenthesis  is  the  same  if  £ and  p are  interchanged.  This  implies 
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A(co)  = \t(ej2Qt'°  '°  + 2 + e-Mk-^)qu 
( =1 


k- 1 k 

+l£  ^L;(2k-f-p)co  + eKp-()(a  + e-Xp-i)a 


+ e-A2k-^-p)(o 


*=1  p=£+l 

In  the  above  expression,  it  is  noted  that  2(k-£)  > 0,  2k-£-p  > 0,  and  p-l  > 0.  The 
expression  also  shows  that  if  ar  is  the  coefficient  of  ^ the  coefficient  of  e'jm  is  also  ar 
for  r > 0.  Therefore,  only  the  nonnegative  powers  will  be  considered.  First  consider  the 
coefficients  of  J®.  Only  the  first  sum  contributes,  and  this  contribution  is  seen  to  be 


l 

4 


*-l 

Yfrw  + <ik.k= 


(=\ 


i 

4 


k- 1 

^2qt(  + 4 qkk 


J= l 


(A4.5) 


the  term  in  brackets  on  the  right  hand  side  of  (A4.5)  corresponds  to  the  first  two  terms  in 
equation  (3.28a).  The  part  of  polynomial  A(co)  containing  positive  powers  of  e 7(0  is  given 
by 


l 

4 


+ 


k-\  k 

^ y/  m-p-oo) 


+ e 


j(p-0  0) 


J=\  p=i+\ 

The  first  term  contains  only  even  powers.  The  second  term  of  the  double  summation  does 
not  contribute  powers  higher  than  k- 1.  Therefore,  consider  the  coefficients  of  e^rw  where 
k<r<2k-2.  If  r is  odd,  only  the  first  term  of  the  double  summation  contributes  anything. 
The  values  of  p and  £ that  satisfy  2 k-p-l  will  contribute  to  this  power.  Therefore, 
p=2k-r-£.  It  is  straightforward  but  involved  to  show  that  the  coefficient  of  this  power  is 


m2 

X2^.2k-r-£ 

e=\ 


where  m2=/c-(r+l)/2.  Therefore, 


ra2 

ar  = £24e,2k-r-(  r odd  (A4.6) 

(=1 

which  is  the  second  formula  of  equation  (3.31). 
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If  r is  even,  but  k<r<2k-2,  then  both  the  first  summation  and  the  first  term  in  the 
double  summation  contribute  to  the  coefficient.  Similar  to  the  case  where  r is  odd,  the 
contribution  of  the  double  summation  is 


m2 

X2^.2k-r-£ 

e=\ 

where  ni2=£-(r)/2.  The  first  sum  will  contribute  when  2 (k-l)  = r,  so  it  follows  that  this 
contribution  is 

9m2  ,m2 


Adding  these  terms  gives 


m2 

— *7m2,/n2  + Z2^',2k-r-£  r even  (A4-7) 

e=i 

which  is  the  first  equation  of  (3.31).  The  coefficients  of  where  I<r<&-1,  given  in 
equation  (3.30a),  are  found  in  a similar  manner.  Using  equations  (A4.5),  (A4.6),  and 
(A4.7)  allows  A(o) ) to  be  expressed  as 


A(co)  = ^ 


k- 1 


2k-2 


4<ik,k  + YMu  + Sar(e7' 

(=1 


JT(a  + e-JTa 


(A4.8) 


r=l 


Defining 


A(^'co)  = 4A(co)  (A4.9) 

results  in  equation  (3.28a).  The  construction  of  the  polynomials  B^7®),  and  D(e7“)  are 
entirely  identical,  and  the  analogues  of  (A4.9)  are 


B(e7(0)  = 4;B(co)  (A4.10) 

D(e7C0)  = -4D(co)  (A4.ll) 

When  equation  (3.26)  is  formed  using  equations  (A4.3)-(A4.4),  and  (A4.9)-(A4. 1 1), 
equation  (3.23)  results.  V 
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5.  Proof  of  Lemma  4.1 


The  constraint  function  is 


<|)(x)  :=  min  Re{a(x,e-/“)} 

we[0,7t]  1 1 

The  function 

RejaCx,^®)} 

with  the  choice  of  a (x,z)  given  by  (4.43)  is  equivalent  to 

Reja(x,e,'a>)j  = l + cT(co)x 

where 

c(to)  = [c](co)  •••  cN(co)  0 •••  0]Te9I2N+1 
The  functions  Cj(co)  are  given  by 

(-aeJ*  + lV 


Cj(o))  = Re 


PJ®  _ a 
V e a j 


for  i = 1,...,N 


Since 


(A5.1) 


(A5.2) 


(A5.3) 


Rejocf^x1  +(1- A.)x2,e7“)|  = 1 + ct(o))[A,x1  +(1-?i)x2 

= X + XcT(co)[x'  j + (1  - X)  + (1  - X)cT((0)[x2] 


= X Re{a(x 1 , <?y“ )}  + (1  - X)  Re{a(x2 , ej(a )} 

for  all  X,  it  follows  immediately  that  (A5.2)  is  an  affine  (and,  therefore,  convex  and 
concave)  function  of  x.  It  is  also  true  that  the  function 

-Reja(x,e;(0)j 

is  affine  (and  thus  convex).  Since  the  maximium  of  a family  of  convex  functions  is  also 
convex,  it  follows  that 


max 

0) 


-Reja(x,e7<0)}J 


is  a convex  function.  Since 
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-max 

CO 


-Re|a(x,e7“)||  = min  Reja(x,e;W)j 

and  the  negative  of  a convex  function  is  concave,  it  follows  that 

<{)(x)  = min  Re|a(x,e-/a))j 

is  a concave  function  of  x.  This  result  also  means  that  all  values  of  x satisfying 


<j)(x)  > y 


is  a convex  set  for  all  real  y. 


The  objective  function  is 

£(x)  = max  ^(x,co) 

CO£[0,7l] 


First  it  is  shown  that  the  function 


(A5.4) 


C(x,m) 


Qp 2 Re{7j  ( eja  )a(x,  eja ) + T2  (eja  )[3(x,  eja )} 
Re{a(x,e710)} 


2 


(A5.5) 


is  quasi-convex,  which  means  that  it  satisfies 

£(X,x'  +(l-^)x2,co)  < max|^(x1,(o),^(x2,co)| 

for  any  X s [0, 1] . This  constraint  is  less  restrictive  than  convexity,  but  it  implies  that  the 
set  of  all  x for  which  £(x,co)<y  is  a convex  set  for  every  value  of  y.  Assume  that 
^(x’,o))  = 8j,  and  £(x2,co)  = 82  for  some  value  of  co,  where  8i  >82  is  assumed 
without  loss  of  generality.  Since  the  set  of  all  x that  satisfy  (|)(x)  > 0 is  convex,  if  it 
assumed  that  x1  and  x2  are  elements  of  this  set,  all  convex  combinations  of  x1  and  x2  are 
also  inside  this  set,  and  are  thus  feasible  points  for  the  objective  function.  Now,  using  the 
affine  nature  of  the  constraint 


C(Ax'  +(l-X)x2,co) 


QfRej 

+(1 

- X)x2,ej(0)  + T2(eJa)$(Xxl  + (1  - X)x2,ej(a)] 

I2 

Re 

a(Xx'  + (l-X)x2,ej(0) 
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Q;/2ARe{r1(^“)a(x\^)  + r2(^“)P(xI,^w)} 

Re{a(Ax'  + (1  - A)x2,e7(0)} 

+ Q;,/2(1-A)R  e{rK^“)a(x2,e^)  + r2(^M)(3(x2,^C0)}||o 


^Q|/2Re 

[ 7j  (*>  )a(x‘ , eja ) + T2  (e7“  )p(  x 1 , ) 

i 

ARe{a(xV“) 

■ + (1  - A) Re{a(x2,e-/0)) 

(1 - A)|qJ/2  Re  | 

[7](^“)a(x2,^“)  + 7'2(^(0)(3(x2 

ARej 

[a(x1,e'“)] 

| + (1  - A)Reja(x2,e7“)j 

1 

A5,  Re| 

a(xV®)] 

> + (1  — A)52  Re | 

[a(x2,e;“)| 

ARejo^x1,^®) 

| + (1  - A)Re{a(x2,e7®)J 

i 

S^ARejafx1,^®) 

} + (l-A)Re{a(xV®)}) 

ARe{ 

a(xV®)} 

+ (l-A)Re{a(xV®)} 

This  implies 


= 5 


l 


^(Ax1 +(1- A)x2,(o)  < max(^(x1,a)),^(x2,co))  = 8] 
and  therefore,  the  function  £( x,co ) is  quasi-convex.  Again,  the  maximum  of  a family  of 
quasi-con  vex  functionals  is  also  quasi-convex,  so  the  objective  function  (A5.4)  is  therefore 
quasi-convex.  V 


6.  Proof  of  Theorem  5.3 

Proof.  Since  from  definition  (5.21)  r(k)  is  the  minimum  of  two  concave  functions,  it 
follows  that  r(k)  must  have  a unique  global  maximum  (Rockafellar  1970).  Hence,  the 
conditions  for  selecting  an  optimizing  controller  in  (5.26a)-(5.26c)  must  yield  a unique 
result.  The  structure  of  the  proof  is  as  follows.  First,  the  conditions  in  (5.26a)  are 
assumed,  and  are  shown  to  lead  to  the  choice  kj  for  the  optimal  controller.  Second,  the 
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conditions  of  (5.26b)  are  assumed,  and  are  shown  to  lead  to  the  choice  of  k2  for  the 
optimal  controller.  Next,  it  is  shown  that  the  conditions  in  (5.26a)  and  (5.26b)  are 
incompatible,  i.e.,  they  cannot  be  simultaneously  satisfied.  Finally,  the  remaining 
compatible  combinations  of  signs  of  a0  + k i (3°  and  a0  + k^P0  are  considered,  along 
with  the  case  where  both  k\  and  k2  are  complex,  and  these  two  possibilities  are  shown  to 
lead  to  the  choice  of  ko  for  the  optimal  controller,  demonstrating  (5.26c)  and  completing  the 
proof. 

Let  Ij  denote  the  range  of  values  of  k where  a0  + kp°  > 0,  and  let  I2  denote  the  range 
of  values  of  k where  a0  + kp°  < 0,  such  that  Ij  u {ko}  u I2  = (—<»,oo).  The  result  (5.26a) 
is  shown  by  demonstrating  that  r(k)  < ri(ki)  on  the  subintervals  Ii  and  I2,  and  at  the  point 
ko-  Assume  that  k]  e and  a0  + kiP°  > 0.  It  is  obvious  that  r(k)  < ri(ki)  for  all  k e I), 
because  r(k)  is  defined  by  the  branch  rj(k)  on  Ij,  and  k]  is  the  maximum  of  n(k).  The 
function  r(k)  = r2(k)  on  I2,  but  from  definitions  (5.22)  and  (5.23)  it  follows  that 
r2(k)  < rj(k)  on  I2,  and  thus  r(k)  < ri(k[)  for  all  k e I2  because  rj (k j ) is  the  maximum 
value  of  rj(k).  Finally,  since  r(ko)  = ri(ko),  it  follows  immediately  that  r(ko)  < ri(kj). 
This  proves  (5.26a).  The  proof  for  (5.26b)  is  analogous.  Next  it  is  shown  that  the 
conditions  a0  + k 1 P°  > 0 and  a0  + k2P°  < 0 cannot  be  simultaneously  satisfied.  The 
definitions  of  k\  and  k2  given  in  (5.24)  and  (5.25)  can  be  rewritten  as  k\  = Ci  - C2P0  and 
k2  = ci  + C2P0,  where  cj,  C2  e 9?,  and  C2  > 0 when  k(  and  k2  are  real.  Thus  if  a0  + k)P° 
= a0  + cip°  - C2(P0)2  > 0,  then  a0  + cip°  + C2(P0)2  >0.  Recognizing  that 
a0  + k2P°  = a0  + C1P0  + c 2( P °) 2 , it  follows  that  if  g°  = a0  + kp°  > 0 at  kj,  then 
g°  > 0 at  k2-  In  a similar  manner,  it  can  be  shown  that  if  g°  < 0 at  k2,  then  g°  < 0 at  ki 
also.  Therefore,  the  only  possibility  not  covered  by  the  conditions  on  the  sign  of  g°  in 
(5.26a)  and  (5.26b)  is  the  case  where  a0  + kiP°  < 0 and  a0  + k2p°  > 0.  To  complete 
the  proof  it  is  shown  that  this  case  leads  to  the  choice  k = ko-  This  pair  of  inequalities 
implies  that  r(kj ) = r2(kj)  and  that  r(k2)  = n(k2).  Therefore,  the  branches  rj (k)  and 
r2(k)  do  not  define  the  curve  r(k)  at  their  maximum  points.  Since  the  function  rj  (k)  has  a 
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unique  maximum,  this  implies  that  the  function  r(k)  does  not  attain  its  maximum  on  Ij, 
since  ^ * 0 everywhere  on  Ij.  Similarly,  r(k)  does  not  attain  its  maximum  on  I2. 
Therefore,  the  only  possible  maximum  is  the  point  ko.  This  implies  that  k = ko.  In  the 
case  that  a0  + ki (3°  = 0,  this  implies  that  the  maximum  of  n(k)  occurs  at  the  point 
kj  = ko.  It  can  be  shown  that  if  g°  = 0 at  ki,  then  g°  > 0 at  k2,  and  thus  the  maximum 
occurs  at  k]  = ko-  Similarly,  if  k2  = ko,  it  can  be  shown  that  the  maximum  occurs  at 
k2  = ko-  In  the  case  where  kj  and  k2  are  complex,  neither  rj(k)  or  r2(k)  has  a real 
maximum.  Again,  r(k)  does  not  acheive  its  maximum  on  Ij  or  I2,  therefore  this  condition 
implies  that  k = ko.  This  proves  (5.26c).  Finally,  from  Corollary  5.1,  it  is  concluded  that 
the  optimal  controller  is  robust  if  r(k)  <0.  V 
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